[q,r]=deconv(b,a) - частное (q) и остаток (r) от деления b на a, так что conv(a,q)+r=b.

residue(b,a) - разложение рациональной функции b(x)/a(x) на элементарные дроби над полем комплексных чисел с выделением целой части. Если a(x) имеет кратные или близкие друг к другу корни, результаты могут быть неверными, поскольку такая задача плохо обусловлена. Плохая обусловленность, т.е. крайне сильная зависимость результата от коэффициентов, иллюстрируется заключительным примером из этой темы.

p=poly(r) - построение полинома по корням, заданным в векторе-столбце r. Для квадратной матрицы r полином p будет ее характеристическим многочленом.

polyval(p,x) - поэлементное вычисление значений полинома p на множестве x, где x может быть как вектором, так и матрицей. Размеры результата совпадают с size(x).

polyder(p) - производная от p.

roots(p) - вектор-столбец, содержащий все корни полинома. Их порядок не определен. Сказанное по поводу неустойчивости результатов для команды residue точно так же относится и к команде roots. Корни полинома вычисляются как собственные значения некоторой матрицы A(p) того же порядка.

Приведем несколько примеров по применению этих команд.

1. Построим график полинома p(x)=x3-x+2 на отрезке -1<=x<=1. Это выглядит так:

p=[1,0,-1,2]; x=-1:.01:1; f=polyval(p,x); plot(x,f), grid

Полином не имеет корней на заданном отрезке. Это подтверждает и команда

roots(p)'

которая даст

ans = -1.5214 0.7607 - 0.8579i 0.7607 + 0.8579i.

2.Разделим предыдущий полином на x-3:

[q,r]=deconv(p,[1,-3])

Тогда

q = 1 3 8, r = 0 0 0 26.

Другими словами, частное q(x)=x2+3x+8, а остаток r=26.

3. Разложим функцию (x-3)/p(x) на элементарные дроби:

1;[r,s,k]=residue([1,-3],p); r', s', k'

Для r' получим вектор из трех компонент r1, r2, r3 :

-0.7607 0.3803 - 0.4289i 0.3803 + 0.4289i,

для s' - также вектор из трех компонент s1, s2, s3 :

-1.5214 0.7607 - 0.8579i 0.7607 + 0.8579i

и k=[] (это означает, что целой части в разложении нет - действительно, у числителя первая, а у знааменателя третья степени). Компоненты векторов r и s означают, что

(x-3)/p(x)=sum(ri/(x-si), i=1:3.

Команда residue работает и в обратную сторону:

2;[q,p]=residue(r,s,k)

восстановит исходные числитель и знаменатель:

q = 0 1 -3 (он получился точно),

p = 1.0000 -0.0000 -1.0000 2.0000 (здесь уже сказались ошибки округления и старший коэффициент не равен 1 автоматически).

4. В заключение приведем сложный пример (Уилкинсон, 1963), показывающий, что иногда, несмотря на хорошую разделенность корней полинома, их вычисленные значения могут очень сильно зависеть от значений некоторых коэффициентов просто потому, что производные корней по этим коэффициентам - очень большие по модулю числа. Такие задачи называются плохо обусловленными и всегда требуют повышенного внимания независимо от того, каким методом они решаются. В то же время они в наибольшей степени стимулируют теоретические исследования по оценке точности машинных вычислений, и пример Уилкинсона - одна из первых классических задач такого рода.

Перейдем теперь к описанию примера. Пусть vn=1:n, где n>1 - целочисленный параметр, и pn=poly(v') - полином с корнями 1:n, которые хорошо отделены друг от друга, а wn=roots(pn) - вектор-столбец с вычисленными корнями полинома pn. Проведем сравнение vn' и wn для различных n. Начнем с n=2:

1;n=2; vn=1:n; pn=poly(vn'); wn=roots(pn); [vn',wn]

и получим ans=1 2 2 1

откуда видно, что элементы в wn нужно упорядочить. Выполняя при n=3 отредактированную строку 1

2;n=3; vn=1:n; pn=poly(vn'); wn=roots(pn); R=[vn',sort(wn)]

найдем R = 1.0000 1.0000

2.0000 2.0000

3.0000 3.0000 .

Появившиеся в первом столбце R цифры 0 как бы "наведены" значениями из второго столбца, и таких мелких шероховатостей у команды format немало. А цифры 0 во втором столбце R говорят о том, что уже появилась погрешность в определении корней. Она, конечно, еще очень мала:

3;((R(:,2)-R(:,1))./R(:,1))'

дает относительную ошибку для корней

ans = 1e-14 *( 0.1110 -0.0444 -0.1184)

- это означает, что в численном результате верны примерно 14 знаков.

Для 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 (это коэффициент при x19). Прибавим к нему 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)=x20+ax19+ ... +20! , a=-210.

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

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

откуда

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

У нас p/a =x19, а полином 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=105 и будет еще больше с ростом x. Поэтому внесение в коэффициент a возмущения 10-7 должно в обязательном порядке заметно сказаться на значениях некоторых корней, каким бы методом они ни находились. Более того, если эти необходимые изменения никак не проявились, метод следует забраковать.

6. Итерации

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

xk+1=F(xk), k=0, 1, 2, ... , x0 задано,

где F - заданное преобразование y=F(x), x0 - как-то выбранное начальное приближение, xk - значение переменной x на k-й итерации, а сама переменная x может быть любой - числом, вектором или матрицей. Предел итераций, если он существует, будем обозначать через X, и тогда уравнение

X=F(X) (1)

должно обязательно выполняться для итерационного преобразования F. Это обстоятельство помогает выбирать различные варианты для F. Все решения X этого уравнения называются неподвижными точками преобразования F(x). Все такие x0, для которых последовательность xk сходится, образуют область сходимости итерационного преобразования F. Cкорость сходимости итераций xk характеризуется поведением числовой последовательности

vk=norm(xk+1-xk)/norm(xk-xk-1),

где норма может выбираться по-разному. Для сходимости xk уже не обязательно, чтобы существовал предел V последовательности vk, но очень часто для сходящихся итераций он существует, и если это так, то пусть a=abs(V). Тогда при a=1 сходимость xk будет крайне медленной, при 0<a<1 это будет сходимость по геометрической прогрессии со знаменателем q=V, а при a=0 последовательность xk будет сходиться быстрее любой геометрической прогрессии. Покажем теперь на простейших примерах, как все это выглядит на деле. Чтобы иметь возможность разнообразия вариантов, задачу возьмем нелинейной. Рвссмотрим уравнение

(x-2)(x-3)=0 или f(x)x2-5x+6=0 c корнями x1=2, x2=3, (2)

построим для нахождения его решений x1=2, x2=3 несколько итерационных преобразований или схем F и проанализируем их работу.

1.Пусть для уравнения (2)

x=F(x), где y=F(x)=(x2+6)/5.

Обязательное условие (1) для преобразования F выполнено, однако при этом переходе появилось дополнительное решение x=inf (). Потеря некоторых или приобретение новых решений часто случается при переходе от исходного уравнения к его итерационной форме. Переходя к нужной нам индексной записи, будем иметь

x(k+1)=F(x(k)), k=1:n,

где начальное приближение x(1) и число итераций n должны быть заданы. Будем накапливать значения x(k) в переменной x, а текущее значение x(k) обозначим через xt. Нужные вычисления реализует строка

1;xt=0; n=100; x=xt; TF='(xt^2+6)/5'; for k=1:n,xt=eval(TF); x=[x,xt]; end, plot(x), grid

Здесь записаны текстовая переменная TF и команда eval (она интерпретирует TF как выражение xt^2+6)/5 и выполняет его). После выполнения строки 1 на графике отобразится 101 значение x(k), включая начальное приближение. Итерации сошлись к px=2. Строка

2;xt=0; n=100; x=xt; TF='xt=(xt^2+6)/5;'; for k=1:n,eval(TF),x=[x,xt]; end, plot(x), grid

произведет те же вычисления, но обратите внимание на то, как в ней записаны TF и eval.

Хотя внешне сходимость x(k) к x1=2 не вызывает сомнений, это в действительности всегда нужно проверять более тщательно. Так как y'=F'(x)=2x/5, то F'(2)=0.8<1, а F'(3)=1.2>1, но итерации, как известно, могут сойтись только к такой неподвижной точке X преобразования F, для которой |F'(X)|1. Отсюда ясно, почему X=2. Однако далеко не всегда можно найти F'(x) аналитически, и поэтому в общем случае приходится использовать вычислительный подход. При известных уже x(k) его реализует строка

3;w=2:n; v=(x(w+1)-x(w))./(x(w)-x(w-1)); plot(v), grid

Страницы: 1, 2, 3, 4, 5