Занятие 6

Материал из theor
Перейти к: навигация, поиск

1. Интегрирование.

1.1 На полупрямой x>=0 вычислить интеграл от функции <math>f(x) = exp(-x^2)</math> и оценить нужный шаг интегрирования dx и верхний предел X.

 dx = 0.1; X=2; x=0:dx:X-dx;
 F=exp(-x.^2)*dx;
 SX=sum(F),plot(x,cumsum(F)),grid;
 dx = 1e-3, X=3, SX = 0.8867;
 dx = 1e-4, X=4, SX = 0.8863;

Точность до единицы третьего знака.

1.2 Теперь вычислим тот же интеграл с помощью ряда и оценим нужное число его членов и X.

Ряд имеет вид суммы по k=0,1,2... членов <math>(-1)^k \frac{X^{2k+1}}{(2*k+1)k!}</math>

 nk = 10; X=2; k=0:nk;
 SX=sum((-1).^k.*X.^(2*k+1)./[1,cumprod(1:nk)]./(2*k+1));

Сначала варьируем nk, затем X (nk = 30, X=3, SX=0.8862).

При каком nk для X=3, и при каком X для nk=30 происходит срыв вычислений? Графики:

 nk = 350; X=3; k=0:nk;
 v=sum((-1).^k.*X.^(2*k+1)./[1,cumprod(1:nk)]./(2*k+1)); plot(v) pause

6zad.png

 nk = 60; X=3; k=0:nk;
 v=(-1).^k.*X.^(2*k+1)./[1,cumprod(1:nk)]./(2*k+1);plot(v) pause

6zad1.png

 nk = 40; X=3; k=0:nk;
 v=(-1).^k.*X.^(2*k+1)./[1,cumprod(1:nk)]./(2*k+1);plot(v) pause

6zad2.png

nk=40 достаточно, для срыва вычислений. Характер потери точности для этих двух случаев существенно различен.

1.3 Вычислить интеграл от функции <math>exp(z-5)</math> по окружности радиуса 2 с центром в точке x=1.

 n=150; fi=linspace(0,2*pi,n);
 z=2*exp(i*fi)+1; plot(z),grid;
 s=sum(exp(z(1:n-1)-5).*diff(z));

Независимость от пути интегрирования.

 n1=60; s1=sum(exp(z(1:n1-1)-5).*diff(z(1:n1)); s2=sum(exp(z(n:-1:n1+1)-5).*diff((n:-1:n1)z)); [s1,s2];

Принцип аргумента:

 N-P=diff(arg(f(z)))/(2*pi).

Для наших функций и кривой <math>f=exp(z-5)</math>

 a=angle(diff(f));
 plot(a);

6zad3.png

 plot(unwrap(a)/2/pi);

6zad4.png

Функция имеет один нуль. Ее график

 plot(f),grid

6zad5.png

2. Логическая задача.

Удаление расположенных по кругу номеров от 1 до n с помощью заданного шаблона s из нулей и единиц. Поясняющий пример :

 n=7, s=[0,1].

Тогда:

 1234567 7246 26 6
 010101  0101 01

Программа запишется так:

 n=7; v=1:n; s=[0,1]; ns=length(s);
 while n>1, r=v(1)*s(1); s=s([2:ns,1]); if r~=0, v(n+1)=v(1); end, v(1)=[]; n=length(v); end, v

Для n=1997 и того же s, v=1946; для n=1997 и s=[1,0,1], v=597.

3. Домашнее задание.

1. Найти нуль функции f из п 1.3 с помощью команды contour.

2. Попробуйте составить алгоритм для задачи из п.2 так, чтобы заметно сократилось время счета.

Нуль функции f для задачи из п.1

 x=-1:0.01:3;
 y=-2:0.01:2;

Преобразуем в массив чтобы можно было вычислять функцию 2-х переменных:

 [X,Y]=meshgrid(x,y);
 v=[0,0]
 f=exp(X+i*Y-5);

Окружность с рад 2 и центром 1

 g=(X-1).^2+Y.^2-4;
 hold on

Contour - дает изолинии z матрицы:

 contour(x,y,f,v)
 pause
 contour(x,y,g,v,'r')
 grid
 hold off
 pause

Итого: 6zad6.png

Алгоритм для задачи из п.2

Попробуем сократить время счета, исключив непрерывное переопределение векторов s и v. Для этого изменим программу, увеличив расход памяти:

 n=1998; s=[1,0,1]'; S=s(:,ones(1,n)); w=find(S==0); s=S(1:w(n-1)); ns=length(s);
 v=[1:n,zeros(1,ns)]; tic, for k=1:ns, r=v(k)*s(k); if r~=0, v(n+1)=v(k); n=n+1; end, end, v(ns+1), toc

Elapsed time is 0.000886 seconds.

Для цикла из п.2

Elapsed time is 0.053445 seconds.

Время счета сокращается в 7-8 раз.