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 интерпретируется системой как полином только тогда, когда он задается в качестве параметра для одной из команд, производящих вычисления с полиномами. Так как в этих командах не проверяется условие an0, надо стараться самим соблюдать его, поскольку иногда это может служить источником ошибок.

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

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.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

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