Ввод векторов

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