и сначала для проверки получим свой старый результат. Но уже при 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. Усовершенствуйте, как можете, обе строки последнего примера и посмотрите, как изменится значение t2/t1.

5. Проверьте, что при вычислении определителя квадратной матрицы M порядка n методом Гаусса нужно выполнить 2n3/3 операций сложения и умножения. Считая их одинаковыми по времени выполнения, найдите, за какое время t выполняется 1 млн операций, взяв n=500 и выдав время счета det(M).

6. Используя предыдущий результат, оцените число арифметических операций для вычисления функций sin и exp.

7. Оцените число операций, которое затрачивается в MATLAB'е для организации одного шага в цикле for-end.

8. Прочтите в help описания команд max, diff, sort и с их помощью попробуйте составить программу, которая для заданной квадратной матрицы M определяла бы, что M является нижней треугольной с ненулевыми элементами на главной диагонали или получается из таковой перестановкой столбцов.

Тема 5. Графический способ решения уравнений

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), при этом сетки по аргументам всегда должны быть прямоугольными.

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