Реферат: Matlab

Для n=10 выполнение строки (ее можно создать из строк 2 и 3)

4;n=10; vn=1:n; pn=poly(vn'); wn=roots(pn); R=[vn',sort(wn)]; R1=(R(:,2)-R(:,1))./R(:,1)

говорит о том, что точность результата постепенно падает. Строка

5;me=max(abs(R1))

дает me=4.2009e-10, т.е. теперь для всех корней верны только 9 знаков (me – max. error). Но корни еще остаются вещественными, поскольку

6; iwn = sum ( abs ( imag ( wn )))

дает iwn=0.

Выполним строку 4 для n=20 и затем строкой 5 найдем максимальную относительную ошибку me=0.0073, так что теперь для всех корней верны только 2 знака, но результат пока еще остается вещественным, поскольку после строки 6 получим iwn=0. Сравним точные и вычисленные корни графически: график

7;plot(R), grid

показывает, что погрешность для некоторых корней уже видна на глаз – для корней 10:17 желтая и фиолетовая линии слегка различаются.

Теперь приступим к самому интересному в данном примере. При n=20 pn(2)= -210 (это коэффициент при x¹⁹ ). Прибавим к нему 1e-7, т.е. внесем в него малое возмущение примерно в 10-м знаке, и повторим расчет с отредактированной строкой 4:

8;n=20; vn=1:n; pn=poly(vn'); pn(2)=pn(2)+1e-7; wn=roots(pn); R=[vn',wn], plot(R), grid

Несмотря на такое малое возмущение в коэффициенте pn(2), некоторые корни стали комплексными. Это видно, во-первых, из выдачи на экране (их мнимые части достигают по модулю 2.7), во-вторых, из того, что теперь строкой 6 получим iwn=18.67, и из графика. Если построить график

9; plot ( R ,'.'), grid

т.е. убрать соединения между точками, результат будет выглядеть более рельефно. На таких графиках режим zoom работает.

Выясним, почему так сильно изменились результаты при внесении столь малого возмущения. Обозначим через p(x) наш невозмущенный полином pn при n=20 и через a его второй коэффициент:

p(x)=prod(x-k),k=1:20, или p(x)=x²⁰ +ax¹⁹ + ... +20! , a=-210.

Тогда для корней x=1:20 по теореме о производной неявной функции будем иметь

¶ p/ ¶ x*dx/da + ¶ p/ ¶ a = 0,

откуда

dx / da = - ¶ p / ¶ a / ¶ p / ¶ x .

У нас ¶ p/ ¶ a =x¹⁹ , а полином ¶ p/ ¶ x находится как polyder(pn) (см. начало этой темы). Поэтому для вычисления dxda на множестве vn наших корней сначала выполним отредактированную строку 4 с n=20

10;n=20; vn=1:n; pn=poly(vn');

а затем строку (на графике значения функции определены только для x=1:20 )

11;dpn=polyder(pn); dxda=-(vn.^19)./polyval(dpn,vn); plot(log10(abs(dxda)),'.'), grid

и увидим, что уже при x=8 dx/da=10⁵ и будет еще больше с ростом x. Поэтому внесение в коэффициент a возмущения 10^-7 должно в обязательном порядке заметно сказаться на значениях некоторых корней, каким бы методом они ни находились. Более того, если эти необходимые изменения никак не проявились, метод следует забраковать.

6. Итерации

Итерации являются с точки зрения программирования одним из самых эффективных способов организации вычислений. Простые итерации в общем случае представляются в виде

x_k+1 =F(x_k ), k=0, 1, 2, ... , x₀ задано,