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.5 d=.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'а:

pi i j eps

Их не следует портить.

Комплексные числа:

q=1+2*i q=1+2i real(q) imag(q) abs(q) conj(q) s=angle(q) (здесь -<s<=).

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:12 toeplitz(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-3 a=1e-5 a=123456789 a=1234567891 a=1+3*i

format long e, 2^.5, format short

Опция format short e позволяет получать ровные столбцы.

Они берутся в кавычки (на букве э на латинском регистре), символ занимает 2 байта. Используются для задания заголовков в числовых выдачах и на графиках, для задания формул и т.д. Можно переводить текстовые переменные в числовые и наоборот. Выполним в командной строке

t='Moscow - столица России' (после дефиса нужно перейти на русский шрифт и затем не забыть снова вернуться на латинский).

Другие типы переменных - ячейки и структуры.

Система help.

help выдает список директорий системы;

help <имя директории> выдает список команд директории;

help <имя команды> выдает описание команды.

type <имя команды> выдает текст команды или программы пользователя, если он составлен в терминах MATLAB'а.

2. Элементы xy-графики

1.Как открывать графическое окно:

figure whitebg zoom on

Теперь построим график функципи y=sin(2x), 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 magenta cyan 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 axis on 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 (при k1/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-inf inf/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)x³y-2xy²+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)=a_nxⁿ+a_n-1x^n-1+...+a₀ задается вектором-строкой p из чисел a_n, a_n-1, ... , a₀,

т.е. коэффициентами, расположенными в порядке убывания показателя степени. Его степень n задавать не надо, поскольку n=length(p)-1; полином может быть и константой - тогда n=0; коэффициенты a_k - любые комплексные числа. Вектор p интерпретируется системой как полином только тогда, когда он задается в качестве параметра для одной из команд, производящих вычисления с полиномами. Так как в этих командах не проверяется условие a_n0, надо стараться самим соблюдать его, поскольку иногда это может служить источником ошибок.

Основные команды для действий с полиномами таковы:

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)=x³-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)=x²+3x+8, а остаток r=26.

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

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

Для r' получим вектор из трех компонент r₁, r₂, r₃ :

-0.7607 0.3803 - 0.4289i 0.3803 + 0.4289i,

для s' - также вектор из трех компонент s₁, s₂, s₃ :

-1.5214 0.7607 - 0.8579i 0.7607 + 0.8579i

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

(x-3)/p(x)=sum(r_i/(x-s_i), 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 (это коэффициент при 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₀ задано,

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

X=F(X) (1)

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

v_k=norm(x_k+1-x_k)/norm(x_k-x_k-1),

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

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

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

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

x=F(x), где y=F(x)=(x²+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) к x₁=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) задерживаются в районе неподвижной точки x₂=3 (они все время уходят от нее, но на графике этого не видно), а затем примерно за 60 итераций монотонно движутся к x₁=2.

(b') xt=3 . Все x(k)=3, а график из строки 3 пуст. Это произошло потому, что при xt=3 F(xt)=15/5=3 получается без ошибок округления.

(b'') xt=3.01 . Пределами для x(k) и v(k) будет inf, так что x₂=3 является неустойчивой неподвижной точкой преобразования F: при малейшем сдвиге x₀ с x₂ в пределе итераций получится либо устойчивая неподвижная точка x₁, либо 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 снова получим x₂ за одну итерацию, а при xt= -3.01 и далее получим inf. Таким образом, при вещественном x₀ итерации сходятся к устойчивой неподвижной точке x₁=2 при -3<x₀<3 и к inf при |x₀|>3. Просчитаем тот случай, когда |x₀|=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=3s*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)), hold on

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 равна 2ⁿ.

2.Теперь представим (2) в виде

x=F(x), где y=F(x)=5-6/x, так что F'(x)=6/x² , F'(3)=2/3<1, |F'(x)|<1 при |x|>2.45.

Следовательно, устойчивая неподвижная точка - это x₀=3. Получим резултат сразу для "всех" вещественных x₀=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, для которой итерации опять идут без ошибок округления. Возникает предположение, что вся комплексная плоскость с выколотой точкой x₁=2 стягивается итерациями в точку x₂=3, хотя не везде |F'(x)|<1. Посмотрим, как это можно увидеть на графике, выполнив программу

2;n=4; fi=-pi:pi/20:pi; xt=4*exp(i*fi); TF='5-6./xt'; z=xt;

3;for k=1:n, xt=eval(TF); z=[z;xt];end, plot(z','.'), axis equal

Начальное приближение (окружность радиуса 4 с центром в нуле) итерируется 4 раза и все результаты выдаются на график. Окружности (на графике их 5) довольно быстро стягиваются в точку x₂=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'')0 и тем самым f'(x')0.

Это и есть итерации по Ньютону для уравнения f(x)=0. При переходе к индексной форме записи для f(x)=x²-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))² .

Так как f'(x_1,2)0, можно использовать эти итерации по Ньютону (часто их называют методом касательной), а поскольку F'(x_1,2)=0, следует ожидать их быстрой сходимости и того, что теперь оба решения x_1,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

Предел итераций при x₀=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 точек.

При x₀=4 предел итераций равен 3, а v(k) успевают подойти к -1 (у v(6) точность уже потеряна). Таким образом, в зависимости от начального приближения x₀=xt итерации могут сходиться к обоим решениям задачи.

Так как теперь преобразование y=F(x) уже не является дробно-линейным, рассмотрим трансформацию комплексной прямоугольной области в процессе итераций. Создадим такую область xt из 41²=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 стягивается итерациями в точку x₁=2, а полуплоскость Re(z)>2.5 - в точку x₂=3: прямая Re(z)=2.5 переходит сама в себя. Точка x₁=2 имеет красный цвет потому, что в нее перешли первые 25 столбцов матрицы xt, а остаток rem(25,7)=4, но 4-й цвет - красный. Далее идут зеленые точки (26-й столбец) и синяя x₂=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 ,

и это есть точки

z₁=2.5-0.5i, z₂=2.5, z₃=2.5+0.5i.

Теперь разберемся, как преобразуется итерациями F прямая Re(z)=2.5. Так как прямая Re(z)=2.5 переходит в себя, ее можно параметризовать с помощью переменной y=Im(z), и пусть y_k=Im(F^k(Re(z)=2.5), k=1, 2,... , т.е. после k итераций y переходит в y_k(y). Ясно, что на k-й итерации в inf перейдут только те точки из всего множества y_k-1, которые равны нулю. Поскольку -inf<y<inf, на 1-й итерации это случится только с одной точкой z₁=2.5 (для нее y=0). Выдадим на график y₁ после 1-й итерации:

10;xt=2.5+(-10:.1:10)'*i; y=imag(xt); n=1; for k=1:n,xt=eval(TF);end, plot(y,imag(xt)), grid

и увидим, что функция y₁(y) пересекает ось абсцисс только два раза (при этом y₁ дважды непрерывно пробегает от -inf до inf), так что на 2-й итерации в inf перейдут только две точки (это уже известные нам z₁=2.5-0.5i, и z₃=2.5+0.5i). Выполним строку 10 с n=2 и снимем с графика строкой

11;q=ginput;q(:,1)'

четыре точки пересечения кривой y₂(y) с осью абсцисс: приближенно это -1.2079 -0.2056 0.2055 1.2079,

тогда как для y₁(y) это были значения -0.5 и 0.5. Каждая переходящая в inf на k-й итерации точка порождает слева и справа от себя две такие точки, которые перейдут в inf на (k+1)-й итерации. Поэтому с ростом k точки y_k(y)=0, с одной стороны, уходят в обе стороны от нуля, а с другой - все чаще появляются в тех местах, где они однажды уже появились. Всего на k-й итерации в inf перейдет 2^k-1 точек. В действительности с ростом k они всюду плотно заполняют прямую Re(z)=2.5, а между ними y_k(y) непрерывно пробегает все значения от -inf до inf. Чтобы быть в этом уверенными, выполним последнюю программу