p(x)=anxn+an-1xn-1+...+a0 задается вектором-строкой p из чисел an, an-1, ... , a0,
т.е. коэффициентами, расположенными в порядке убывания показателя степени. Его степень n задавать не надо, поскольку n=length(p)-1; полином может быть и константой – тогда n=0; коэффициенты ak – любые комплексные числа. Вектор p интерпретируется системой как полином только тогда, когда он задается в качестве параметра для одной из команд, производящих вычисления с полиномами. Так как в этих командах не проверяется условие an¹0, надо стараться самим соблюдать его, поскольку иногда это может служить источником ошибок.
Основные команды для действий с полиномами таковы:
conv(p,q) – произведение полиномов p и q. Название команды происходит от слова convolution (свертка), поскольку коэффициенты произведения действительно получаются как компоненты свертки векторов p и q.
[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.8579i0.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 (µ). Потеря некоторых или приобретение новых решений часто случается при переходе от исходного уравнения к его итерационной форме. Переходя к нужной нам индексной записи, будем иметь