Реферат: Matlab
Министерство образования РеспубликиБеларусь
Учреждение образования
«Гомельский государственныйуниверситет им. Ф. Скорины»
Математический факультет
Кафедра МПУ
MATLAB
Реферат
Исполнитель:
Студентка группы М-53
Гумарева Л.С.
Гомель 2004
Введение
MATLAB – матричнаялаборатория – наиболее развитая система программирования для научно-техническихрасчетов, дополненная к настоящему времени несколькими десятками более частныхприложений, относящихся к вычислительной математике, обралботке информации,конструированию электронных приборов, экономике и ряду других разделовприкладной науки. Изучение MATLAB'а по фирменной документации, которая теперьприлагается на инсталляционном компакт-диске, занимает у начинающихпользователей слишком много времени не только из-за необходимости читать ее наанглийском языке со специфическим слэнгом, но, главным образом, ввиду неизбежногодля таких руководств последовательного и достаточно формального изложениябольшого объема информации, а имеющиеся на русском языке пособия в основномследуют этому стереотипу. Даже для опытного специалиста по расчетам накомпьютере такое изучение сопряжено с неоправданно большими затратами труда.
MATLAB предназначенпрежде всего для программирования численных алгоритмов. Он разрабатывается ужеболее 15 лет и возник на основе более ранних прикладных пакетов LINPACK иEIGPACK, созданных в 1970-е гг. в США, и в свою очередь повлиял на появлениетаких систем, как MathCad, MAPLE и Mathematica. Совершенствование системыMATLAB происходило как в связи с достижениями в вычислительной математике, таки в связи с изменениями в архитектуре персональных компьютеров и развитиемобщесистемных средств. Со временем MATLAB был дополнен целым рядом ужеупоминавшихся приложений (toolboxes), далеко раздвинувших границы егоприменимости. Далее речь пойдет лишь о ядре MATLAB'а, которое мы будем называтьсистемой, и конкретно о ее версии 5.2, выпущенной фирмой MathWorks в январе1998 г.
MATLAB – системапрограммирования высокого уровня, работающая как интерпретатор и включающаябольшой набор инструкций (команд) для выполнения самых разнообразныхвычислений, задания структур данных и графического представления информации.Команды эти разбиты на тематические группы, расположенные в различныхдиректориях системы. Теперь в системе около 800 команд, и примерно половина изних вполне доступна начинающему пользователю. Команды с большим возможнымобъемом вычислений написаны на С, но много и таких команд, которые представленыв терминах этих первых. Поэтому система оказывается почти открытой дляпользователя. Имеются большие возможности для вывода двумерной и трехмернойграфики и средства управления ею. Пользователь может без особых затрудненийдобавлять свои команды и писать программы в терминах уже существующих команд;несколько сложнее делать это в рамках Фортрана и С. Можно обмениваться даннымис программами на этих языках, а из них обращаться к системе. Краткость инаглядность программирования и исключительные возможности визуализациирезультатов делают систему очень эффективной при поисках и апробации новыхалгоритмов, при проведении разовых расчетов и в учебном процессе, поскольку ееможно осваивать без предварительного знакомства с основами программирования ивыполнять такие сложные примеры, которые невозможно делать с использованиемдругих систем.
Документация по системе иее приложениям содержит много тысяч страниц, и поэтому естественно встаетвопрос о том, как ее осваивать. Работа с системой требует определеннойматематической подготовки, так что обучение можно начинать на втором курсевуза. Основные сведения о системе изложены в руководствах /1/ – /2/: /1/ – этоучебник с описанием вычислительных возможностей и архитектуры системы, /2/ –описание ее графических возможностей. Конечно, можно читать подряд /1/, /2/ ипри необходимости обращаться за уточнениями к команде help или справочнику /3/,в котором описаны почти все команды. Но гораздо более эффективным, на нашвзгляд, является изложение основных вычислительных процедур с помощью наиболееупотребительных команд системы. Именно так мы и познакомимся с MATLAB'ом, аточнее примерно с 30-40 его командами. После этих занятий пользоватьсядокументацией /1/ и /2/ будет гораздо легче.
Приложений к системе мыздесь не касаемся, а изучать их можно только после предварительногоознакомления с ней, а также с тем разделом науки, которому посвящено конкретноеприложение. Отметим только, что большинство приложений означают дляпользователя просто расширение списка доступных ему команд. Очень удобно то,что вся документация по системе и приложениям находится на компакт-диске, скоторого происходит их установка, и при желании она может быть размещена такжеи на винчестере.
Для работы с системойдостаточно иметь компьютер PC 486 с оперативной памятью хотя бы 16 Mb и сустановленными на нем системами Windows 95 и MATLAB 5.2. В действительностиMATLAB может работать и с друогими операционными системами, такими, например,как Macintosh, Unix и OS/2.
За рубежом вышло ужедостаточно много учебных пособий по системе, но на русский язык ни одно из нихпока не переводилось и даже в центральных библиотеках их теперь нет из-засокращения финансирования. Изданные у нас пособия (например, /4/ – /12/) восновном следуют руководствам /1/ – /3/, тогда как нам представляется полезнымдать менее формальное введение в предмет, опираясь прежде всего на интуициюслушателя.
1. Переменные
Переменные могут быть числовыми, текстовыми идругих типов. У нас будут только числовые (это во всех деталях) и текстовые(совсем немного). Название переменной начинается с латинской буквы, далее могутбыть буквы и цифры (не более 31 символа). Строчные и прописные буквы здесьразличаются.
1. Числовые переменные.Это числа, векторы, матрицы и многомерные массивы. В компьютере все числапредставлены примерно с 16 десятичными знаками, под каждое вещественное числоотводится 8 байтов, под комплексное – 16.
1.1 Ввод чисел
Целые числа. В системеони не выделяются явно. Наберем и выполним отдельно каждую команду:
a=2 a=2.0 a=2;a=1:6 b=1:20 c=10:-2:5
Командное окно. Команднаястрока. Редактирование командной строки. Буфер исполненных команд. Как выбиратьинформацию из командного окна и из буфера исполненных командных строк. Нельзядопускать совпадения имени переменной с именем какой-либо команды.
Вещественные числа.Выполним по отдельности следующие команды:
d=0.5:0.3:2.5d=.5:.3:2.5 d=.5+1:.3-.1:2.5*2 length(d)
d(end) d(end-2)d(1) d(0) d(2:7) d(7:-1:2) d(150)
f=linspace(1.5,30,143); length(f)
Индексы всегда начинаютсясо значения 1. Команды набираются на малом латинском регистре. Возможнамногопараметричность команд.
Диапазон вещественныхчисел:
realmax realmin
Другие константыMATLAB'а:
piijeps
Их не следует портить.
Комплексные числа:
q=1+2*i q=1+2ireal(q) imag(q) abs(q) conj(q) s=angle(q) (здесь -p<s<=p).
q=1+2*i;r=3;fi=0:.01:pi; z=q+r*exp(i*fi); plot(z) Это верхняя полуокружность.
1.2 Ввод векторов
Векторы-строки:
a=1:6 linspace(1,6,10)
Векторы-столбцы:
a=(1:6)' linspace(1,6,10)'
Операторы .' и ' :
y1=linspace(1,6,4)';y2=y1; y=y1+i*y2; y.' y'
Команды linspace и: применимы для задания тольковещественных векторов.
Ввод матриц. A(i,j) — элемент из i-й строки и j-го столбца. A(k) – k-й элемент таблицы,вытянутой в столбец.
A=[1,2;3,4]A=[1;2,3;4] A(2,2) A(3) A(5) size(A) A(3,4)=10 size(A)
A(5)=6 size(A)A(22)=3 A=A(:) A(22)=3 size(A) [m,n]=size(A)
A=reshape(1:24,4,6)size(A) A([1,end],:)=[] A(:,[1,end])=[] size(A)
1.3Некоторыеспециальныематрицы
m=3;n=4; eye(m,n)eye(m) eye(n) ones(m,n) ones(m) ones(n) zeros(m,n)
rand(m,n) rand(m,n) rand('state',0)rand(m,n) rand(m) Эторавномерное распределение на интервале (0, 1).
randn(m,n) randn('state',100) Это нормальное распределение, у негомат.ожидание=0, дисперсия=1
v1=1:4 v2=7:12toeplitz(v1,v2) toeplitz(v1)
1.4 Некоторыепростыекоманды
A=reshape(1:24,4,6)triu(A) triu(A,0) triu(A,2) triu(A,-1) tril(A)
v=1:5 diag(v)diag(v,2) diag(v,-1)
diag(A) diag(A,2)diag(A,-1)
A=reshape(1:24,4,6)rot90(A) rot90(A,2)
Выдачи на экран. Командаformat с различными опциями.
В обычном формате(forrmat short) выдается 5 знаков, для целых чисел 9 знаков, порядки изменяютсяот -308 до +308. В полном формате(format long e) 16 знаков.
a=2 a=.001 a=1e-3a=1e-5 a=123456789 a=1234567891 a=1+3*i
format longe, 2^.5, format short
Опция format short e позволяет получать ровные столбцы.
Они берутся в кавычки (набукве э на латинском регистре), символ занимает 2 байта. Используются длязадания заголовков в числовых выдачах и на графиках, для задания формул и т.д.Можно переводить текстовые переменные в числовые и наоборот. Выполним вкомандной строке
t='Moscow — столицаРоссии' (последефиса нужно перейти на русский шрифт и затем не забыть снова вернуться налатинский).
Другие типы переменных –ячейки и структуры.
Система help.
help выдает список директорий системы;
help<имя директории> выдает список команд директории;
help<имя команды> выдает описание команды.
type<имя команды> выдает текст команды или программыпользователя, если он составлен в терминах MATLAB'а.
2. Элементы xy-графики
1.Как открыватьграфическое окно:
figurewhitebgzoomon
Теперь построим графикфункципи y=sin(2px), 0<=x<=5, выполнивстроку
x=0:1e-3:5;y=sin(2*pi*x); plot(y) plot(x,y) ,grid
Использование режима zoom:
k=100; y=sin(2*pi*k*x); plot(y)
2.Автоматическоечередование цветов. Теперь будем, как правило, нумеровать строки.
1;x=linspace(0,1,20);k=.1:.1:.8; y=k'*x; plot(x,y)
Здесь определяетсявектор-строка x=0:20, затем вектор-строка k из 8 угловых коэффициентов, далееполучается матрица y=k'*x какпроизведение вектора-столбца k' на вектор-строку x. Строки этой матрицы состоятиз точек соответствующих прямолинейных отрезков. Наконец, строятся графики этихотрезков как функций от x – первая нижняя линия (она желтая) соответствуетk=.1, последняя, тоже желтая, – для k=.8. Мы видим, что цвета, которых всего 7,чередуются циклически в таком порядке (под русскими английские названия):
желтый фиолетовый голубойкрасный зеленый синий белый
yellow magentacyan red green blue white
Вызовем строку 1 иотредактируем в ней команду plot:
1;x=linspace(0,1,20);k=.1:.1:.8; y=k'*x; plot(x,y,'g.')
т.е. добавим там третий(текстовой, ибо он в апострофах) аргумент. Все кривые на рисунке станутзелеными (green), а линии будут изображаться отдельными точками. Аналогичноупотребляются и другие цвета из этого списка – по первой букве. В текстовомаргументе может быть до трех символов. Для изображения точек графика помимо. употребляютсяеще: — -. * x o + и некоторые другие символы.
3.Графики в полярныхкоординатах:
x=1:.01:3;nx=length(x); r=x.^2; fi=linspace(0,5*pi,nx); polar(fi,r)
4.Еще один пример – легкостроятся многозначные функции:
x=0:.1:6*pi;y=cos(x); plot(x,y) plot(y,x)
5.Управление осями:
axis off axison axis ([-10,10,-5,20]) axis auto axis equal axis square
Размеры осей можнозадавать и для трехмерной графики, но цвета в ней используются дляхарактеристики величины ординаты и команда zoom там не работает.
3. Простые примеры,иллюстрирующие эффективность MATLAB'а
1. Суммирование.Найдем при заданном n частичную сумму ряда s(n) = 1/k^2, k=1:n. Для этоговыполним строку
1;n=100;k=1:n; f=k.^(-2); plot(cumsum(f)), [sum(f),pi^2/6] =1000
Команда cumsum(f)подсчитывает все частичные суммы s(k) от f(1:k) для каждого k от 1 до n, такчто на графике можно наблюдать процесс формирования нужной нам величины. Вконце строки выдается численный и точный результаты:
ans = 1.6350 1.6449.
Полагая n=1000, получим
ans = 1.6439 1.6449,
т.е. ошибку в 1 единицу4-й значащей цифры.
Сходимость не всегдастоль очевидна, как на этом графике. Чтобы в этом убедиться, усложним нашпример: при заданных m>1 и n найдем частичную сумму ряда s(m,n) =sum(1/k^m), k=1:n (при m=1 получается уже расходящийся гармоническийряд). Для проведения вычислений отредактируем строку 1:
2;m=2; n=1000; k=1:n;f=k.^(-m); plot(cumsum(f)), sum(f)
=1.5 =1e4
=1.2
и сначала дляпроверки получим свой старый результат. Но уже при m=1.5 у нас, глядя награфик, нет полной уверенности в достижении сходимости. Это тем более так приm=1.2: для n=1000 ans=4.3358, а для n=1e4 ans=4.7991. Факт сходимости ряда приm=1.01 нельзя установить численно из-за низкой скорости его сходимости.
Чтобы лучше запомнитьдействие команды cumsum, вычислим ò(x/sin(x))dx, xÎ[0, 3]. Подинтегральная функция f=x/sin(x) не имеет в нулеособенности, и поэтому достаточно выполнить строку
3;n=100;h=3/n; x=h/2:h:3-h/2; f=x./sin(x); plot(h*cumsum(f)), grid, sum(h*f) =1000
т.е. аппроксимировать f всерединах интервалов (эти точки x называют полуцелыми в отличие от концовсчетных интервалов – целых точек). Сравнение ответа ans = 8.4495 играфика наводит на мысль о том, что пока сходимость еще не достигнута, но приn=1e3 ans = 8.4552, так что при n=1e2 со сходимостью в действительности все впорядке, а возрастание функции h*cumsum(f) на правом конце происходит из-зароста там функции f – это можно увидеть, выполнив
4;plot(f)
Для матрицы A команды sumи cumsum работают вдоль столбцов (значит, по первому индексу), а для вектора –вдоль него независимо от того, строка это или столбец. Чтобы провестисуммирование для матрицы A вдоль ее строк, нужно выполнить sum(A,2), т.е.указать для выполнения команды второй индекс. Это правило относится ко многимкомандам MATLAB'a и к многомерным матрицам тоже – по умолчанию имеется в видупервый индекс, а в противном случае нужно всегда указывать, по какому индексудолжна работать команда, и это указание не сохраняется для последующих команд.
2. Произведения.Аналогично суммированию с помощью команд prod и cumprod вычисляются иобрабатываются произведения. Например, найдем Õ(1-1/k^2), k=2:100 (при k®¥ Õ®1/2), выполнив строку
1;n=100; k2=(2:n).^2; a=1-1./k2; cp=cumprod(a); cp(end), plot(cp/.5), grid
Результат cp(end) =0.5050 говорит о том, что сходимость здесь не очень быстрая. Это видно и изграфика, на котором представлена относительная ошибка результата. Обратитевнимание на названия переменных k2=k^2 и cp=cumprod(..): при выборе именпеременных всегда нужно стремиться к тому, чтобы эти имена хоть как-то отражалисуть дела (это особенно важно при написании больших программ, где многопеременных).
При вычислениипроизведений можно выйти за числовую шкалу. Найдем, например, для каких k можнонайти k!.. Ясно, что максимально допустимое km вряд ли больше 200, так чтострока
2;n=200;k=1:n; kf=cumprod(k); plot(kf)
должна дать ответ на нашвопрос. Из-за быстрого возрастания kf и ограниченной разрешимости дисплея (этоне более 0.5% от максимального значения на графике) мы видим всего одну точкуkf(km), перед которой, как нам ошибочно кажется, идут нули и за которой идутчисла inf (infinity), вообще никак не представленные нарисунке. Точно так же графика обходится и с переменной NaN (not a number), и это обстоятельство может быть иногда полезным.Переменная NaN возникает в таких ситуациях:
0/0 inf-infinf/inf
Переменные inf и NaN (они получаются со знаком) можно использовать впрограммах. Для определения km выполним строку
3;sum(isinf(kf))
в которой isinf(kf)выдаст 1 на тех позициях вектора размеров kf, где элементы kf есть inf, и 0 наостальных позициях. Поскольку ans=30, km=n-30=170, что можно было бы получить исразу, выполнив строку
4;km=sum(isfinite(kf))
где isfinite отмечает теэлементы числовой переменной, которые отличны от inf и NaN. При выходепроизведения за числовую шкалу для сомножителей можно использовать команды
log (взятие натуральногологарифма),
log10 (взятие десятичногологарифма),
abs (взятие модуля),
sign (взятие знака,выдающее 1, 0 и -1).
3. Логические задачи. Обычно при освоении программированиялогические действия даются труднее арифметических. Приведем здесь два простыхпримера задач логического характера.
1. Напишем строку длянахождения общих элементов двух векторов:
x=1:20;y=15:30; [X,Y]=meshgrid(x,y); v=X(X==Y)
2. Второй пример несколько сложнее, иначинающие изучать MATLAB обычно пытаются решить его с помощью циклов for-end,что совершенно неправильно. Взяв на сторонах единичного квадрата по 200интервалов, определим, сколько точек получившейся таким образом сетки попадаетвнутрь вписанной в него окружности. Нужная программа имеет вид
1;tic,x=0:1/200:1; [X,Y]=meshgrid(x); M=abs(X+i*Y-.5-i*.5)<1/2; s=sum(M(:)),t1=toc
и даст ответ s=31397точек, t1=0.16 сек, тогда как строка для циклов for-end
2;tic,s=0;w=1:201; for I=w,for J=w,if norm([x(I),x(J)]-.5)<.5,s=s+1; end,end,end,s ,t2=toc
дает то же самое s иt2=7.47 сек, так что t2/t1=46. Это лишний раз говорит о том, что нужно разумноподходить к использованию операторов языка программирования.
4. Графический способрешения уравнений
1. Простой пример: найтикорни уравнения x*sin(x^2)=0 на отрезке [0,3]. Программа:
1;x=0:.01:3; f=x.*sin(x.^2); plot(x,[f;0*f]), grid
2;ginput
В команде ginput точкаснимается нажатием левой клавиши мыши, Enter – выход из ginput.
Проверим это другимспособом:
3;nx=length(x); w=1:nx-1; x(find(f(w).*f(w+1)<0|f(w)==0)) Отв: 0, 1.77, 2.5.
Эту строку можноупростить:
4;nx=length(x); w=1:nx-1; x(f(w).*f(w+1)<0|f(w)==0)
Матрицы и векторы сэлементами 0-1.
2. Сложный пример –неявные функции. Построим график неявной функции f(x,y)ºx3y-2xy2+y-0.2=0,x,y=[0, 1]. Это выполнит программа
1;h=.02; x=0:h:1; [X,Y]=meshgrid(x); f=X.^3.*Y-2*X.*Y.^2+Y-.2;
2;v=[0,0];contour(x,x,f,v), grid
На графике зеленая линия(справа она двузначная) представляет искомый результат. Область в первомквадранте между этими кривыми обозначим через G. Эту задачу совсем непростосделать в других системах программирования прежде всего потому, что вычислениеобразующих линии уровней точек – в общем случае очень сложная процедура.
Выясним, какой знак имеетf в области G, для чего выполним
3;mesh(x,x,f.*(f>0))
Это пример трехмерной,т.е. xyz-графики. В ней цвет используется дляизображения амплитуды (значения z),
изменяясь с ростом zот темносинего через голубой,зеленый и желтый до темнокрасного.
Вычислим площадь S этойобласти:
4;S=h^2*sum(f(:)>=0) (S=0.7296).
Для h=0.01 выполнимстроку 1, затем строку 4 и получим S=0.7204, а для h=0.005 найдем S=0.7152. Приинтегрировании всегда естественно делать такие проверки.
Выясним, какой объемзаключен между поверхностью f(x,y) и областью G, где f(x,y)>=0. Для этогоснова возьмем в строке 1 h=0.02 и вычислим
5;V=h^2*sum(f(f>=0)) (V=0.1268)
Для h=0.01 V=0.1235, адля h=0.005 V=0.1219. Теперь не нужно писать f(:), поскольку f(f>=0) естьвектор.
Конечно, эти результатыприближенные (с точностью до 1 — 2%), но отметьте, как быстро и просто они былиполучены. Такие приемы можно применять для решения достаточно широкого кругазадач.
Выполним строку
6;C=contour(x,x,f); clabel(C)
которая зашлет числовуюинформацию о графике в матрицу C ипостроит график, выбрав значения уровней автоматически. Из матрицы C можно последовательно выбирать всекривые.
Обобщения. Графическим способом можно решатьсистемы уравнений и уравнения в комплексной плоскости. Команда contour3 строит линии уровней для функций f(x,y,z), при этом сетки по аргументамвсегда должны быть прямоугольными.
5. Полиномы
По степени применимости,по разнообразию и качеству соответствующих команд скалярные полиномы –следующие за матрицами математические объекты в MATLAB'е. Полином
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.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.00002.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=8dx/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'; fork=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;'; fork=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
Из графика видно, чтоv(k) действительно сходятся к a=0.8,как и должно быть согласно теории, которая здесь выглядит очевидной.
Теперь будем варьироватьначальное приближение, выполняя строки 1 и 3.
(a) xt=1.5, 1.9, 1.99, 1.9999.При последних двух значениях xt на графике из строки 3 появятся осцилляциисправа, так как здесь разности x(k+1)-x(k) и тем самым значения v(k) теряютточность: x(k) уже сошлись со многими знаками.
(a') xt=2. График изстроки 2 вообще пуст, потому что тогда все v(k)=0/0=NaN.
(b) xt=2.01, 2.5, 2.9, 2.99,2.9999. В последнем случае x(k) вначале довольно долго (до k=20) задерживаютсяв районе неподвижной точки x2=3 (они все время уходят от нее, но награфике этого не видно), а затем примерно за 60 итераций монотонно движутся к x1=2.
(b') xt=3. Все x(k)=3, аграфик из строки 3 пуст. Это произошло потому, что при xt=3 F(xt)=15/5=3получается без ошибок округления.
(b'') xt=3.01. Пределамидля x(k) и v(k) будет inf, так что x2=3 является неустойчивойнеподвижной точкой преобразования F: при малейшем сдвиге x0с x2в пределе итераций получится либо устойчивая неподвижная точка x1,либо inf – приобретенная неподвижная точка. Проварьируем этот сдвиг:
(с) xt=3+1e-15, 3+1e-14, 3+1e-8. В 1-м варианте ухода xt не происходит, во 2-м тенденция ухода уже зародиласьи он неизбежно произойдет с увеличением числа итераций, в 3-м уход проявился вполной мере уже на 100 итерациях.
(с') xt=-2.5, -2.9, -2.99,-3, -3.01, -100, 100. При xt=-3 снова получим x2 за одну итерацию,а при xt= -3.01 и далее получим inf. Таким образом, при вещественном x0итерации сходятся к устойчивой неподвижной точке x1=2 при -3<x0<3и к inf при |x0|>3. Просчитаем тот случай, когда |x0|=3.Этот расчет выполняется строкой (редактируем строку 1)
4;n=100;fi=-pi:pi/20:pi; xt=3*exp(i*fi); TF='(xt.^2+6)/5'; for k=1:n,xt=eval(TF);end,plot(xt,'.')
На графике (онкомплекснозначный) видны 4 точки: точка z=3, точки z=3±s*i с малым s>0 и точка z=2. Чтобы разобраться в результате, выполнимстроку 4 с n=1000 и, сделав окно MATLAB'а полным, выдадим
5;imag(xt')
Образы первой и последнейточки начальной окружности слегка отличаются от z=2, а ее 21-я точка (онасоответствует fi=0) есть z=3. Это получилось потому, что sin(pi) и sin(-pi) неравны нулю в точности. Снова сделаем окно MATLAB'а небольшим и выполним строку4 с радиусом 3.01, а затем строки
6;sum(isnan(xt))
7;find(isnan(xt))
и получим, что точкипервоначальной xt с номерами 1, 21, 41 обращаются в inf (здесь nan=inf/inf).Для радиуса 5 их будет 21, для радиуса 5.6 их будет 41, т.е. все.
Границу областисходимости обычно трудно исследовать аналитически, даже в таком простомпримере, как этот, и сама она не определяется из условия |F'(x)|=1 или |x|=2.5.Условие |F'(X)|<1 гарантирует устойчивостьнеподвижной точки X преобразованияF(x), условие |F'(X)|>1 являетсядостаточным для ее неустойчивости, а в случае |F'(X)|=1 она может быть как устойчивой, так и неустойчивой.Неустойчивые неподвижные точки сравнительно редки в вычислительных задачах, ноих полезно иметь в виду при исследовании численных алгоритмов. В нашем случаемы установили только, что множество |x|£3 принадлежит области сходимости итераций F, но вовсе неописали границу этой области. Мы установили также, что при |x|³5.6 итерации сходятся к inf, и поэтому inf является устойчивойнеподвижной точкой F. Чтобы получить представление о границе областисходимости, выполним строки
8;r=3:.1:5.6;z=[];n=100; for kr=r,xt=kr*exp(i*fi); for k=1:n,xt=eval(TF); end,z=[z;xt]; end
9;zn=isnan(z);Z=r'*exp(i*fi); plot(Z(zn)), axis equal
и увидим приближенныйграфик границы – это вовсе не окружность. Комагнда axis equal выбирает по осям одинаковыймасштаб (это действует только на текущий plot; у команды axis много другихопций). Чтобы построить границу аккуратнее, выполним строку 4 с шагом по fi, равным pi/180, затем строку 8 с r=3:.02:5.6 и строку 9. Получим график границы, выполнивстроки
10;zn=isnan([z;z(end,:)]); zn=diff(zn)~=0; plot(Z(zn)), holdon
11;y=3*exp(i*fi);plot(y,'m'), axis equal, hold off
и увидим отличие областисходимости от круга |z|£3.
Найдем те точки, которыеперейдут в z=3 на первых 10 итерациях. Для этогопридется рассмотреть обратное к Fпреобразование x=±(5y-6)^.5 и сделать с ним 10 итераций, задав начальное значение y=3. Строка
12;n=10;yt=3; for k=1:n,yt=(5*yt-6).^.5; yt=[yt,-yt];end, plot(yt,'.')
показывает, что при этомполучится практически та же линия, что и в строке 10. Удивительно, что этоверно и для других точек z изобласти сходимости преобразования F, но с некоторыми вкраплениями внутрь области сходимости. Брать n большим нельзя, так как в концедлина вектора yt равна 2n.
2.Теперь представим (2) ввиде
x=F(x), гдеy=F(x)=5-6/x, так что F'(x)=6/x2, F'(3)=2/3<1, |F'(x)|<1 при|x|>2.45.
Следовательно, устойчиваянеподвижная точка – это x0=3. Получим резултат сразу для «всех»вещественных x0=xt:
1;n=100;xt=-5:.5:5; x=xt; TF='5-6./xt'; for k=1:n,xt=eval(TF);end, plot(x,xt,'.')
Все пределы равны 3,выпадает только неустойчивая точка xt=2, для которой итерации опять идут безошибок округления. Возникает предположение, что вся комплексная плоскость свыколотой точкой x1=2 стягивается итерациями в точку x2=3,хотя не везде |F'(x)|<1. Посмотрим, как это можно увидеть на графике,выполнив программу
2;n=4; fi=-pi:pi/20:pi; xt=4*exp(i*fi); TF='5-6./xt'; z=xt;
3;fork=1:n, xt=eval(TF); z=[z;xt];end, plot(z','.'), axis equal
Начальное приближение(окружность радиуса 4 с центром в нуле) итерируется 4 раза и все результатывыдаются на график. Окружности (на графике их 5) довольно быстро стягиваются вточку x2=3. Применим zoom, чтобы посмотреть малые окружности.
Сделаем разрезы наначальной окружности:
4;n=4;fi=-pi:pi/20:pi*.75; fi(end-4)=[]; xt=4*exp(i*fi); TF='5-6./xt'; z=xt;
и снова выполним строку3. Мы увидим, что 1-я итерация меняет направление обхода окружности напротивоположное, остальные – нет. Выполним строку 4 с n=20 и затем строку 3. Спомощью zoom дойдем до самой малой окружности (она белого цвета) и убедимся,что она лежит правее точки z=3 примерно в полосе от 3.0001 до 3.0004.
Потерь и приобретенийдругих неподвижных точек здесь нет.
3. Решим нашу задачуметодом Ньютона. Если x'' удовлетворяет уравнению f(x'')=0, а x' находитсявблизи x'' и используется для приближенной аппроксимации f(x''), то
0=f(x'')=f(x')+f'(x')(x''-x') или x''=x'-f(x')/f'(x') при условии, что f'(x'')¹и тем самым f'(x')¹0.
Это и есть итерации поНьютону для уравнения f(x)=0. При переходе к индексной форме записи для f(x)=x2-5x+6и f'(x)=2x-5 из (2) получим
x(k+1)=x(k)-f(x(k))/f'(x(k))или y=F(x), где F(x)=x-f(x)/(f'(x), F'(x)=2f(x)/(f'(x))2.
Так как f'(x1,2)¹0, можно использовать эти итерации поНьютону (часто их называют методом касательной), а поскольку F'(x1,2)=0,следует ожидать их быстрой сходимости и того, что теперь оба решения x1,2будут устойчивыми неподвижными точками итерационного преобразования F.Неподвижной точкой будет и inf, но она «не наша», и природа ее, какмы увидим ниже, гораздо сложнее.
Выделим TF в отдельнуюстроку и проведем наши стандартные вычисления
1;TF='xt-(xt.^2-5*xt+6)./(2*xt-5)';
2;xt=0;n=100; x=xt; for k=1:n,xt=eval(TF); x=[x,xt];end, plot(x), grid
3;w=2:n;v=(x(w+1)-x(w))./(x(w)-x(w-1)); plot(v), grid
Предел итераций при x0=0равен 2, а сходимость имеет порядок выше первого, поскольку v(k) до потери имиточности успели подоойти к нулю. Чтобы определить порядок сходимости,отредактируем строку 2 на предмет оценки квадратичной сходимости:
4;w=2:n;v=(x(w+1)-x(w))./(x(w)-x(w-1)).^2;plot(v(1:6)), grid
Теперь до потери точностиv(k) успевают подойти к 1 (у v(7) точность уже потеряна), так что сходимостьитераций здесь квадратичная. Напомним, что точность теряется у v(k), но не уx(k). Чтобы увидеть потерю точности у v(k), выполним строку 4, выдавая в plot 7 точек.
При x0=4предел итераций равен 3, а v(k) успевают подойти к -1 (у v(6) точность ужепотеряна). Таким образом, в зависимости от начального приближения x0=xtитерации могут сходиться к обоим решениям задачи.
Так как теперьпреобразование y=F(x) уже не является дробно-линейным, рассмотрим трансформациюкомплексной прямоугольной области в процессе итераций. Создадим такую областьxt из 412=1681 комплексных точек и проитерируем ее (при этом вкомандном окне будет много сообщений о делении на нуль):
5;x=-10:.5:10;[X,Y]=meshgrid(x); xt=X+i*Y; n=100; for k=1:n,xt=eval(TF);end, plot(xt,'.')
На графике будут обарешения задачи и некоторое множество точек на прямой Re(z)=2.5. Сделаем только10 итераций и вставим выдачу графика на каждой итерации:
6;x=-10:.5:10;[X,Y]=meshgrid(x); xt=X+i*Y;n=10; for k=1:n,xt=eval(TF); plot(xt,'.'),pause(0), end
Конечный результат тотже, что и при 100 итерациях, но видна динамика его формирования.
ПолуплоскостьRe(z)<2.5 стягивается итерациями в точку x1=2, а полуплоскостьRe(z)>2.5 — в точку x2=3: прямая Re(z)=2.5 переходит сама в себя.Точка x1=2 имеет красный цвет потому, что в нее перешли первые 25столбцов матрицы xt, а остаток rem(25,7)=4, но 4-й цвет – красный. Далее идутзеленые точки (26-й столбец) и синяя x2=3, поскольку rem(41,7)=6 и6-й цвет – синий. Но всегда лучше проверить такие выводы численно: найдем
7;z=xt(:); [sum(abs(z-2)<.1), sum(abs(z-3)<.1), sum(abs(real(z)-2.5)<.1)]
(=1025=41*25), (=615=41*15),(=38<41 на 3).
Так как потерь вматрице-таблице xt MATLAB не допускает, выясним, во что перешли эти 3 точки — вinf или NaN:
8;z=xt(:,26);[sum(isinf(z)), sum(isnan(z))] =0, =3 .
Индексы этих трех точекиз 26-го столбца находятся командой
9;find(isnan(z))' =20 21 22 ,
и это есть точки
z1=2.5-0.5i, z2=2.5, z3=2.5+0.5i.
Теперь разберемся, какпреобразуется итерациями F прямая Re(z)=2.5. Так как прямая Re(z)=2.5 переходитв себя, ее можно параметризовать с помощью переменной y=Im(z), и пусть yk=Im(Fk(Re(z)=2.5),k=1, 2,…, т.е. после k итераций y переходит в yk(y). Ясно, что наk-й итерации в inf перейдут только те точки из всего множества yk-1,которые равны нулю. Поскольку -inf<y<inf, на 1-й итерации это случитсятолько с одной точкой z1=2.5 (для нее y=0). Выдадим на график y1после 1-й итерации:
10;xt=2.5+(-10:.1:10)'*i; y=imag(xt); n=1; fork=1:n,xt=eval(TF);end, plot(y,imag(xt)), grid
и увидим, что функция y1(y)пересекает ось абсцисс только два раза (при этом y1 дваждынепрерывно пробегает от -inf до inf), так что на 2-й итерации в inf перейдуттолько две точки (это уже известные нам z1=2.5-0.5i, и z3=2.5+0.5i).Выполним строку 10 с n=2 и снимем с графика строкой
11;q=ginput;q(:,1)'
четыре точки пересечениякривой y2(y) с осью абсцисс: приближенно это -1.2079 -0.2056 0.2055 1.2079,
тогда как для y1(y)это были значения -0.5 и 0.5. Каждая переходящая в inf на k-й итерации точкапорождает слева и справа от себя две такие точки, которые перейдут в inf на(k+1)-й итерации. Поэтому с ростом k точки yk(y)=0, с одной стороны,уходят в обе стороны от нуля, а с другой — все чаще появляются в тех местах,где они однажды уже появились. Всего на k-й итерации в inf перейдет 2k-1точек. В действительности с ростом k они всюду плотно заполняют прямуюRe(z)=2.5, а между ними yk(y) непрерывно пробегает все значения от-inf до inf. Чтобы быть в этом уверенными, выполним последнюю программу
12;hy=1.0033e-4; xt=2.5+(-1:hy:1)'*i; w=2:length(xt); n=100; fork=1:n,xt=eval(TF); end
13;pzn=sum(sign(xt(w-1)).*sign(xt(w))<0)
которая зафиксирует 674перемены знака (pzn) у функции y100(y) для -1£y£1 с hy=1.0033e-4 (из-за недостаточной малости hy далеко невсе они учтены). Для -10£y£-8 с hy=1.0033e-4 pzn=675. Длясимметричного относительно точки y=0 отрезка 8£y£10 с тем же hy получим pzn=702. Ошибки при последовательномвычислении Fk(z) на прямой Re(z)=2.5 будут ничтожными, так чтозначения yk(y) вычисляются с высокой точностью, но от шага hy числонайденных перемен знака в yk(y), конечно, зависит: при k=100 числовсех нулей функции yk(y) равно 2100-1=6.3383e29.
Сделаем краткое резюмеотносительно неподвижной точки inf преобразования F. Ее следует рассматриватькак неустойчивую, поскольку любая ее окрестность с разрезом Re(z)=2.5 вкомплексной плоскости с ростом k «уходит» от нее, а из этого разрезавсе больше точек попадает в нее, и эти ее дискретные прообразы все плотнее (впределе всюду плотно) заполняют этот разрез, но все их множество, поскольку оносчетно, имеет меру нуль. Отсюда следует, что на прямой Re(z)=2.5 почти всюду несуществует теоретического предела итераций. Из-за ошибок округления, хотя они ималы, прообразы inf почти никогда непопадают в inf при вычислениях и потому ведут себятак же (т.е. хаотически), как и не прообразы.
Заметим теперь, чтопорядок действий в F можно переставить:
F(x)=x/2+1.25+0.25/(2x-5)и F(x)=(x2-6)/(2x-5)– это две другие математически эквивалентные записи F.Выполним строку
14;TF='xt/2+1.25+.25./(2*xt-5)';
а затем строки 5 и 7 иполучим тот же результат, что и выше. Но строка
15;TF='(xt.^2-6)./(2*xt-5)';
и строки 5, 8 дадуттолько устойчивые неподвижные точки преобразования F и, как и раньше, три точкиNaN, так что с такой TF мы не увидели бы рассмотренной выше картины. Выполнимстроку 6 с n=100 и увидим, как все происходит доконца.
Пример 3 показывает, что,пользуясь довольно простыми командами MATLAB'а и руководствуясь здравымсмыслом, можно проанализировать достаточно тонкие математические детали задачи.
7. Системы линейныхалгебраических уравнений
В MATLAB'e имеетсянесколько команд для решения таких задач. Рассмотрим здесь наиболееупотребительную из них – оператор \ (обратный слэш или деление слева). Длячисел a\b=a-1b в отличие от того, что a/b=ab-1. Дляневырожденной квадратной матрицы A и вектора-столбца f вектор-стобец u=A\f естьрешение системы линейных уравнений Au=f, т.е. в записи по аналогии с числамиu=A-1f, но матрица A-1 в действительности не вычисляется,а система Au=f решается методом исключения Гаусса (это значительно дешевле), скоторым каждый из нас как-то знаком еще со школы. Правая часть f может быть иматрицей, и тогда каждый столбец матрицы u=A\f есть решение системы для соответствующегостолбца из f, т.е. командой A\f можно сразу решить много систем, и поступатьтак всегда выгодно, ибо тогда матрица A будет разлагаться на множители пометоду Гаусса только один раз. Если матрица A не квадратная, вектор u=A\f естьрешение системы Au=f в смысле метода наименьших квадратов, но мы не будем здесьразбирать этот вариант оператора \, поскольку не должны сколько-нибудь заметновдаваться в детали численных методов. Таким оброазом, оператор \ – одна изсамых мощных команд системы. А поскольку решение систем Au=f с квадратнымиматрицами – наиболее часто встречающаяся в вычислительной практике задача,необходимо иметь представление хотя бы об одном из методов ее реализации.
1.Решим и проанализируемна точность систему Au=f, выполнив следующую строку:
1;x=1:.1:5;m=length(x); A=toeplitz(exp(x)); ut=sin(x)'; f=A*ut; u=A\f;
Здесь задаетсявектор-строка x=1:.1:5 из 41 элемента, затем по вектору exp(x) командойtoeplitz вычисляется квадратная матрица A размеров 41*41, далее в видевектора-столбца задается точное решение ut, по нему строится правая часть f и,наконец, находится численное решение системы Au=f. Выясним вид A, построивграфики
2;mesh(A) plot(A(1,:))plot(A(20,:)) plot(A(41,:))
Теперь выполним команды
3;plot(ut) plot(f) plot([ut,u])
и увидим, что f сильноотличается от u, а точное и численное решения графически совпадают(изображающая u фиолетовая линия накрыла желтую линию ut). Иногда нужносравнивать сильно разномасштабные кривые (у нас это u и f). Для такогосравнения следует провести нормировку каждой из них на свой абсолютныймаксимум. В нашем случае график
4;plot([u/max(abs(u)),f/max(abs(f))])
позволяет изобразитьдинамику изменения кривых u и f относительно друг друга.
2.Команда det вычисляетопределитель. Найдем
1;max(abs(u-ut)) (=9.7167e-13) det(A) (=4.1058e-9)
откуда видно, что u и utсовпадают примерно в 12 знаках, хотя определитель det(A) системы на первыйвзгляд очень мал. Тем не менее попробуем решить нашу систему по правилуКрамера, обозначив такое решение через uc:
2;d=det(A);uc=zeros(m,1); for k=1:m,uc(k)=det([A(:,1:k-1),f,A(:,k+1:m)])/d; end,plot([ut,uc])
Здесь при k=1 A(:,1:k-1)=[],а при k=m A(:,k+1:m), поскольку 1:0=[] и m+1:m=[], так что все матрицы Ck=[A(:,1:k-1),f,A(:,k+1:m)],k=1:m, сформированы правильно. Из графика видно, что так найденное uc сновасовпадает с точным решением ut. Теперь
3;max(abs(uc-ut)) (=9.6934e-13) т.е. погрешность практически не изменилась. Длительное времясчитали, что прималых d систему решать численно нельзя.
3.В действительноститрудности численного решения системы связаны с возможной близостью к нулюсобственных значений матрицы A, которые определяются как все решения l1, ..., lmхарактеристического уравнения det(A-lE)=0, где E=eye(m) – единичная (с нулями вне главнойдиагонали) матрица порядка m. За последние 10 лет были созданы достаточнохорошие программы для нахождения собственных значений. В MATLAB'е это делаеткоманда eig, которую мы и применим сейчас:
1;sp=eig(A);plot(sp), [y,yi]=min(abs(sp)), [z,zi]=max(abs(sp)), c=z/y
Из графика видно, что всесобственные значения вещественны (по оси абсцисс отложены значения индекса) иникак не упорядочены. Ближе всего к нулю (=0.136) третье из чисел lk, дальше всего (=898) – 41-е, а модуль их отношенияc(A)ºc=6.6040e+3 называется числомобусловленности (condition number) матрицы A и отражает гораздо болееглубокие ее свойства, чем величина ее определителя (которая по существу вообщеничего не отражает ни в практическом, ни в теоретическом плане). При обращении[V,D]=eig(A) в диагональной матрице D расположены lk (раньше они были в векторе-столбце sp), а в матрице Vв виде столбцов – соответствующие им собственные векторы vk сединичной среднеквадратичной нормой
sum(vk2))1/2=1, k=1:m.
Команда eig (и целый рядоснованных на ней) имеет неплохую точность при решении спектральных задачсредней сложности и является очень информативной. В нашем примере естьотносительно близкие друг к другу lk.Чтобы в этом убедиться, выполним строку
2;ssp=sort(sp);v=1:m-1; di=diff(ssp); min(abs([di./ssp(v);di./ssp(v+1)]))
(нормировка разностейделается на каждый член пары) и получим 0.0044, т.е. есть такая пара, у которойсовпадает более двух знаков, так что наш пример в спектральном плане не совсемтривиальный. Теперь сделаем матрицу A вырожденной путем дублирования ее строк ипосмотрим, как это отразится на результатах eig:
3;r=zeros(1,m-1);for k=1:m-1, v=[1:k,k,k+2:m]; C=A(v,:); r(k)=min(abs(eig(C))); end, plot(r)
Результаты следуетпризнать неплохими – max(r)=1e-13 и есть даже несколько чистых нулей.
Чтобы получитьпредставление о собственных векторах преобразования A, выполним строку
4;[V,D]=eig(A);D=V'*V; D(1:m+1:m^2)=0; mcv=max(abs(D(:)))
и получим mcv=7.5460e-16.Это означает, что собственные векторы с высокой степенью точности ортогональны,как и должно быть для симметричной матрицы A.
4.Теоретическая оценкапогрешности при решении линейных систем. Объясним смысл числа обусловленностиc(A), который имеет фундаментальное значение для вычислительных методовлинейной алгебры. Пусть
A*u=f¹0, A*du=df¹0, re(f)=max(abs(df))/max(abs(f)).
Тогда u¹0 и du¹0 ввиду невырожденности A. Будемрассматривать df как ошибку в f (ввиду линейности задачи это не приводит кпотере общности – считать вектор df малым не нужно). Оценим величину
re(u)=max(abs(du))/max(abs(u))для всех допустимых f и df,
т.е. возникающую прирешении нашей системы относительную ошибку, которая является весьма общейхарактеристикой матрицы A. Если y, z и c – числа, определенные в строке
1;sp=eig(A);plot(sp), [y,yi]=min(abs(sp)), [z,zi]=max(abs(sp)), c=z/y
(это строка 1 изпредыдущего примера), то
max(abs(du))£z-1max(abs(df)),
поскольку du=A-1uи z-1 – максимальное по модулю собственное значение для A-1,причем равенство для max(abs(du)) теоретически возможно. Аналогично
max(abs(u))³y-1max(abs(f)),
поскольку y-1 — минимальное по модулю собственное значение для A-1, и знакравенства здесь также возможен. Объединяя оценки для числителя и знаменателяre(u), будем иметь
re(u) £z-1/y-1(max(abs(df))/max(abs(f)))или re(u) £(y/z) re(f), т.е. re(u) £c*re(f),
причем равенствообязательно достигается хотя бы для одной пары f и df. Другими словами, числообусловленности c(A)есть точная верхняя граница роста относительной ошибкипри решении нашей системы: если относительная ошибка правой части равна re(f),то относительная ошибка решения re(u) не превзойдет ее более чем в c(A) раз.
Ввиду важности этогорезультата полезно понимать его и на пальцах: df переходит в du с возрастаниемкомпонент не более чем в z-1 раз, а f переходит в u с уменьшениемкомпонент не более чем в y-1 раз, так что
re(u) £(z<sup/>-1/y-1) re(f).
Мы установили смысл числаc(A), исходя из максимум-нормы векторов f, когда norm(f)=max(abs(f)), но болееточное рассмотрение показывает, что это верно и для среднеквадратичной нормы
norm(f)=(sum(abs(f.^2))^0.5.
Команда cond(A) вычисляетзначение c(A) для квадратной невырожденной матрицы A, но в cond спектр spищется не для A, а для B=A'A, ибо тогда вся спектральная задача для Bстановится заведомо вещественной и, более того, все ее собственные векторывзаимно ортогональны (последнее обстоятельство существенно уменьшаетсреднеквадратичную ошибку при этих сложных вычислениях). Для получения наших yи z нужно лишь извлечь квадратный корень из y(B) и z(B) – это также делается вкоманде cond(A). Для нашего примера c(A)=
2;cond(A) (=6.6040e+3)
и совпадает со значением cиз строки 1, поскольку A — симметричная матрица.
Несмотря на внешнююпростоту, понятие числа обусловленности было четко сформулировано только всередине 1960-х гг., т.е. примерно через 15 лет после появления первых ЭВМ, ашироко использоваться стало еще позже. Для вычислительных методов оно оказалосьважнее понятия линейной зависимости, которое теперь некоторым образомвыражается через него, но мы уже не будем здесь этого уточнять.
5.Практическая оценкапогрешности. Число обусловленности может быть не всегда подходящим для оценкипогрешности. Это так при большом n – тогда решение спектральной задачи вкоманде cond будет слишком дорогим. Это так и в том случае, когда ошибка dfспецифически неравномерна и имеет характерный профиль – тогда желательно иметьи профиль для поточечных ошибок du(k), чего, конечно, никак не может датьcond(A). В таких случаях приходится прибегать к непосредственному моделированиюошибок, которое состоит в задании некоторого множества случайных возмущенийправой части и решении всех таких систем с последующим поточечным описаниемграниц получившихся решений. Проиллюстрируем этот способ оценки ошибки на нашемпримере.
Восстановим этот пример:
1;hx=.1;x=1:hx:5; m=length(x); A=toeplitz(exp(x)); ut=sin(x)'; f=A*ut; u=A\f;
(здесь введен шаг hx).Пусть ошибка
df(x) в f(x) естьравномерно распределенная случайная величина,
по модулю непревосходящая значения g(x)=1e-4*òf(t)dt в пределах от 1 до x.
В строке
2;g=1e-4*hx*abs(cumsum(f));rand('state',0); n=30; V=(1:m)'*ones(1,n); plot([f/max(abs(f)),g/max(g)]), grid
задается профиль g=g(x)максимально допустимой ошибки в точке x, приводится в исходное состояниесчетчик случайных чисел rand, задается число n=30 возмущений правой части f иматрица V размеров m*n из продублированного n раз столбца (1:m)'. Графикпоказывает, что самые большие ошибки там, где abs(f) мало. В строке
3;F=f(V)+g(V).*(2*rand(m,n)-1);U=A\F;plot(U), pause, plot([g(V).*(2*rand(m,n)-1),g])
создается матрица Fвозмущенных правых частей из m строк и n столбцов путем прибавления кразмноженному n раз вектору f возмущений в нужных границах, решаются все системыAU=F и выводятся на график все решения U и все возмущения вместе с их границей.Из последнего графика видно, что возмущения заданы правильно. В строке
4;ua=max(U,[],2); ui=min(U,[],2); plot([ui,ut,ua])
путем обработки U вдольстрок находятся поточечные границы решений (ua – верхняя, ui – нижняя) истроится график точного решения и этих границ.
Результат кажется намплохим потому, что на графике
5;plot(F)
возмущений просто невидно (масштаб f забил их), а в U они отразились слишком сильно. Но формальноон не противоречит оценке, связанной с числом обусловленности c(A).Действительно, без учета профиля df максимальная ошибка me вычисляется как
6;sp=eig(A); me=min(abs(sp))^(-1)*max(g)
и равна 0.6064, тогда как
7;max([ua-ut;ut-ui]) (=0.5180)
и лишь немного подрастетс увеличением n (при n=100 это 0.5540), т.е. не превзойдет значения me, что исвидетельствует о формальной согласованности обоих методов оценки погрешности.Чтобы окончательно избавиться от ощущения какой-то несогласованности нашихметодов, применим критерий cond(A) для среднеквадратичной нормы. Для этогог намнужно пройтись по всем столбцам k=1:n (сейчас n=30) матриц U и F и вычислить
max((norm(du)/norm(u)) / (norm(df)/norm(f)) ).
Запишем это в виде(максимально используйте карман при наборе строк)
8;max((sum((U-ut(V)).^2).^.5./sum(U.^2).^.5)./(sum((F-f(V)).^2).^.5./sum(F.^2).^.5))
и получим 3.6484e+3, чтоменьше c(A)=6.6040e+3, так что и здесь мы не обнаружили противоречия.
Если подойти к оценкепогрешности упрощенно, построив график
9;gm=max(g);plot(A\[f-gm,f,f+gm])
то все три линии на немпрактически совпадут. Таким способом можно моделировать ошибку, еслипреобразование A-1 монотонно, т.е. если при f1<f2обязательно u1<u2 или, наоборот, обязательно u1>u2.Но у нас это не так: на графике (эта строка получается из строки 9)
10;gm=100*max(g);plot(A\[f-gm,f,f+gm])
желтая линия (онасоответствует решению для правой части f-gm<f) то выше, то ниже фиолетовой.Именно из-за немонотонности преобразования A-1 получается такойзаметный разброс в U. С помощью этого приема можно быстро выяснить монотонностьпреобразования y=Q(x) (не обязательно линейного), что далеко не всегда удаетсяопределить теоретически: если все отклики Y=Q(X), x-a<X<x+a, a>0,лежат между Q(x-a) и Q(x+a), то преобразование Q монотонно, и нужно лишь взятьзначение a таким, чтобы результат был виден на графике (для этого нам пришлосьувеличить max(g) в 100 раз).
6.Посмотрим, какизменятся результаты нашего примера при увеличении m — порядка матрицы A.Выполним, не меняя смысла задачи, отредактированную строку 1 предыдущегопримера
1;clearall, hx=.01; x=1:hx:5; A=toeplitz(exp(x)); ut=sin(x)'; f=A*ut;u=A\f;plot([ut,u]), c=cond(A)
Здесь шаг hx уменьшен в10 раз, так что теперь порядок m=401 – довольно высокий; c(A)= 6.0804e5 возрослопочти в 100 раз, т.е. обусловленность A заметно ухудшилась (примерно в 102раз), но
2;max(abs(u-ut)) (=1.3153e-9)
еще достаточно мал, хотяи возрос примерно в 103 раз, т.е. больше, чем c(A). Такоерасхождение с теорией как бы предупреждает о том, что даже при сохранениисмысла задачи увеличение ее размерности не позволяет автоматически применятькритерий числа обусловленности к оценке ошибок округления. К выбору числа mнужно всегда относиться с повышенным вниманием.
Чтобы получитьпредставление о собственных векторах преобразования A, выполним строку
3;[V,D]=eig(A);D=V'*V; m=length(x); D(1:m+1:m^2)=0; mcv=max(abs(D(:)))
и получим mcv=2.2985e-15,т.е. степень ортогональности остается удивительно высокой. Жордановы клеткипорядка выше первого могут быть тогда, когда mcv(A)>0.99.
Мы рассмотрели этотпример так подробно, чтобы показать исключительно высокие возможности MATLAB'а втом, что касается анализа результатов.
Заключение
MATLAB – высокоуровневаясистема программирования, позволяющая резко сократить затраты труда припроверке алгоритмов и проведении прикидочных расчетов. Возможность проведениябольших расчетов на MATLAB'е определяется в основном теми затратами времени, накоторые может пойти пользователь: здесь приходится выбирать между легкостью инаглядностью программирования и представления результатов, с одной стороны, изатратами времени на счет – с другой. Система очень удобна для освоения иапробации численных методов, что мы и хотим показать здесь прежде всего. Именнопоэтому она рекомендуется как одна из основных для физиков и многих другихестественно-научных специальностей в ведущих американских университетах.Детальное освоение любой большой программной системы – это достаточнодлительный процесс, основу которого составляют индивидуальная работа, и нашизанятия призваны дать лишь первоначальный импульс этому процессу в отношенииMATLAB'а. Темы 2 – 4 представляют сравнительно элементарное введение, а востальных рассматриваются более сложные примеры, показывающие, как можноиспользовать программные и графические возможности системы для исследованиячисленных алгоритмов.
Литература
1. Using MATLAB. Version5.2. The Mathworks, Inc., 1997. 531 p. MATLAB 5.2 Product Family New Features.Version 5.2. The Mathworks, Inc., 1998. 202 p.
2. Using MATLAB Graphics.Version 5.2. The Mathworks, Inc., 1997. 372 p.
3. MATLABFunctions Reference (Volumes 1 and 2).Version 5. The Mathworks, Inc., 1998. 819 p., 586 p.
4. Дьяконов В.П. Справочник поприменению системы PC MatLab. М., Физматлит, 1993 112 с.
5. Потемкин В.Г. Система MATLAB.Справочное пособие. М., «Диалог-МИФИ», 1997. 350 с.
6. Гультяев А. MATLAB 5.2. Имитационное моделирование всреде Windows. СПб, «Коронс-принт»,1999, 288 с.
7. Дьяконов В.П., Абраменкова К.В. MATLAB 5. Система символьной математики.М., Нолидж, 1999, 633 с.
8. Лазарев Ю.Ф. MATLAB 5.х. Киев, Изд. группа BHV, 2000, 384 с. («Б-кастудента»).
9. Медведев В.С., Потёмкин В.Г. Control System Toolbox.MATLAB 5 для студентов.М., «Диалог-МИФИ», 1997, 287 с.
10. Потёмкин, В.Г. Введение в MATLAB. М., «Диалог-МИФИ», 2000,350 с.
11. Потёмкин, В.Г. Система инженерныхрасчетов MATLAB 5.х. В 2-х томах. М.,«Диалог-МИФИ», 1999, 366 с., 304 с.
12. Рудаков П.И., Сафонов В.И.Обработка сигналов и изображений. MATLAB 5x. М., Диалог-МИФИ", 2000, 413 с.(«Пакеты прикладных программ»).