Примечание: Укажите ошибку вычислений при выполнения п.8.
Как ее исправить? См. п.10.
9.
Построить в декартовых координатах кривую, заданную параметрически на [
;] (рисунок 8): t
cos v
t sin u
)
t
(
z
)
t
(
y
)
t
(
x
)
t
(
r
, «тильда»– первые три компонента одноименных векторов.
>>
MyCurve1=@(t) u(1)*sin(t)+v(1)*cos(t);
MyCurve2=@(t) u(2)*sin(t)+v(2)*cos(t);
>> % Обратите внимание, что мы здесь не делаем матриц, используя облегченный вызов функции- рисовалки – поэтому не можем применять индексацию по осям
>> MyCurve3=@(t) u(3)*sin(t)+v(3)*cos(t); ezplot3(MyCurve1,MyCurve2,MyCurve3,[-pi,pi])
26
Рисунок 8 – Координатная кривая
10. Создать киноленту, отображающую эволюцию поверхности треугольников при изменении третьей стороны треугольника.
>> % Формируем сетку и область переменной времени, она же равна с
>> a=0:0.01:1;b=0:0.01:1;c=0:0.01:2;[X,Y]=meshgrid(a,
b);
>> % Последовательно создаем графические объекты: рисунок, пространство осей, поверхность
(пока пробную)
>>
HFig1=figure(1); axes();
Hsur=surf(X,Y,ones(length(X))); set(HFig1,
‘PaperType’, ‘A4’, 'Color','y');
>> % Модифицируем параметры объектов. Убираем сетку, делаем цветовую гамму неподвижной (важен порядок вызова установки пределов третьей оси и цветовой легенды
>> colormap('jet'); set(gca,'Box','on'); caxis([0 0.5]); colorbar; set(Hsur,'LineStyle','none')
27
>> % Основной цикл – смена данных построения.
На лету формируем массив из скриншотов-картинок
>> for j=1:length(c) set(Hsur,'Zdata',geron1(X,Y,c(j))); MyMovie(j) = getframe(HFig1); end;
>> % Проигрываем кино 5 раз подряд. Записываем avi-файл, его можно проиграть стандартными средствами Windows.
>> movie(MyMovie,5); movie2avi(MyMovie,'mymovie1.avi','fps',4)
Задание:
1.
Сформировать случайную комплекснозначную матрицу А размера 3×3, вычислить матричную экспоненту
A
e
B
. Решить систему линейных уравнений
T
)
1
,
1
,
1
(
Az
2.
Сформировать, не используя циклов, матрицу вида
ij ij b
a
A
B
AB
и выделить из нее вещественную (D1) и комплексную
(D2) часть
3.
Вычислить с помощью анонимной функции электрический потенциал системы точечных электрических зарядов, заданных координатами (D1,D2), в точке с координатами (0,0). Заряды положительные, по модулю равны 1.
4.
Построить контурный график потенциала на площади
]
2
;
2
[
]
2
;
2
[
. Выполнить цветовую легенду в логарифмическом масштабе. В графическом редакторе оформить график красиво.
5.
Построить в полярных координатах
2D-график напряженностей поля (velocity plot).
6.
Построить 3D-график п.4.
7.
Создать анимацию контурного графика в п.4. при движении координат зарядов по закону:
))
jt i
sin(
a
1
(
2
D
2
D
))
it j
cos(
a
1
(
1
D
1
D
ij ij ij ij
,
0Примечание: поскольку графика требует времени для своего
28 построения, то использовать команду pause для задания задержек.
Более предпочтительный путь – команда drawnow до взятия фрейма.
29
Лабораторная работа №3
Проведение вычислений с помощью М-файлов
1. Создание и отладка М-файлов
MATLAB является интерпретирующим языком непосредственных вычислений, т.е. выражения, которые вы вводите, интерпретируются и вычисляются. Поскольку все вычисления в
MATLAB выполняются с двойной точностью, формат вывода может управляться с помощью следующих команд (рисунок 9).
Рисунок 9 – Команды управления форматом вывода
После вызова одного из приведенных выше форматов он сохраняется до вызова другого. Команда format compact подавляет большинство пустых строк, позволяя большее количество информации вывести на экран или страницу. Она не зависит от других команд формата. Команда eval чрезвычайна полезна, позволяя обмениваться между процедурами (особенно, внешними) строкой- параметром, которая принуждает систему к выполнению заданной этой строкой команды (или серии команд). Напомним, что MATLAB различает большие и маленькие буквы в именах команд, функций и переменных, а спецзнаки могут быть записаны через обратный слэш.
Для простых операций удобен интерактивный режим, но если вычисления нужно многократно выполнять или необходимо реализовывать сложные алгоритмы, то следует использовать m- файлы MATLAB (расширение файла состоит из одной буквы m).
Познакомимся со script-m-файлами (или сценариями)– текстовыми файлами, содержащими инструкции на языке
MATLAB, подлежащими исполнению в автоматическом пакетном режиме.
Создать такой файл удобнее с помощью редактора системы MATLAB
(не забудем способ создания, исходя из окна истории команд). Он
30 вызывается из командного окна системы MATLAB командой меню
File/New/M-file (или самой левой кнопкой на полосе инструментов, на которой изображен чистый белый лист бумаги). Записанные в script-файлы
команды будут выполнены, если в командной строке ввести имя script-файла
(без расширения).
Переменные, определяемые в командном окне и переменные, определяемые в сценариях, составляют единое рабочее пространство системы
MATLAB, причем переменные, определяемые в сценариях, являются глобальными, их значения заместят значения таких же переменных, которые были использованы до вызова данного script-файла. Текст реально имеющихся m-файлов (системы MATLAB, например polar, или ваших собственных) можно просмотреть с помощью команды type <имя_функции>.
Однако, большинство m-файлов являются файлами-функциями, т.е. программами. При том же порядке вызова в отличие от простого линейного скрипта, хотя бы тот и содержал, например, циклические конструкции, нормальный m-файл содержит тело, окаймленное следующими декларациями: function [out1, out2, ...] = funname(in1, in2,
...)
End
Тело может быть представлено в таком формате:
=
%Строки основного комментария, доступного по команде help
% Далее три четыре основных повторяющихся типа конструкций
% Объявление собственных переменных и команды
MATLAB
<переменная=выражение>; <команда>
% Объявление подфункций (пишется обычно с
31 отступом) function [out]=nested_funname (in)
<операторы;> end;
% Конструкции языка типа for, if,….
If <условие>
<операторы;> end
В заголовке явном виде пишутся выходные аргументы, что несколько отличается, например, от Си-программ (там слово function равнозначно procedure). Сходные декларации имеют место для подфункций и пишутся внутри тела, которое по сути тождественно пакету команд. MATLAB рассматривает m-файл как одну «главную»
(primary) функцию, причем ее имя обязано совпадать именем m- файла; внутри этой функции могут располагаться декларации локальных (вложенных, nested) функций, доступных только изнутри данного m-файла. Ниже главной функции в редких случаях могут располагаться подфункции, но доступ к ним возможен только из главной – поэтому их лучше все-таки реализовывать как вложенные.
Внутри m-файлов можно ссылаться на другие m-файлы, в том числе и на самого себя рекурсивно. Тем самым реализуются, хотя и не полностью, принципы процедурного программирования. Переменные в функциях являются по умолчанию локальными, но в версиях 4.0 и выше разрешено объявлять требуемые переменные глобальными
(global). Знак “;” получает двойное толкование: с одной стороны – это подавление вывода в консоль, с другой – привычный разделитель между операторами языка.
Как и в библиотечных функциях MATLAB, используется принцип избыточного задания аргументов. Это позволяет варьировать действия в зависимости от формы и количества аргументов (см. nargin, nargout и особенно параметр декларации функции varargin). Аналогично, доступ к результатам выполнения может быть частичным. При вызове функции для нее создается стек памяти, в который помещаются помимо локальных переменных копии одноименных переменных вызова. Для простого интерфейса, использующего только консоль, предусмотрены операторы disp,
32 input, error. Как и в большинстве языков/сред программирования,
MATLAB обладает стандартными операторами ветвления, выбора и цикла (кроме goto). Особенность for заключается в использовании массива/матрицы в качестве индекса, при этом
счетчик последовательно пробегает значения, быть может вещественные, элементов/столбцов матрицы. Например, если вам необходимо выполнить <оператор> только для тех элементов матрицы, которые больше 3, то удобно это сделать следующим образом: for i=find(A>3) <оператор> end;
Редактор/отладчик предоставляет как средства редактирования текста m-файла, так и средства пошаговой его отладки. Один из способов вызова редактора – вызов из командной строки MATLAB с помощью команды edit. Редактор, используемый в системе, имеет синтаксическую раскраску, т.е. слово или символ по мере ввода приобретают тот цвет, который соответствует их типу. С помощью пункта меню Tools-Fonts можно настроить такие важные параметры, как используемый шрифт. Это особенно важно для работы с русским текстом, поскольку не все шрифты правильно воспроизводят русский текст. Редактор имеет стандартный набор возможностей (запуск M- файла, расстановка точек останова – breakpoints), так и специальные: переключение в Cell Mode (режим ячеек), публикация html, пункт
Evaluate Selection, который позволяет вычислять значение выделенного выражения и помещать результат в консоль, пункт Edit-
Paste to Workspace. В режиме ячеек легко отлаживать и читать программу; с той же целью введена пиктограмма Show Functions, позволяющая перескакивать по заголовкам вложенных функций.
Более детальные сведения можно прочесть, нажав пункт меню редактора Help – Using the M-file Editor.
Справа от редактируемого файла находится линейка номеров строк, а слева (рядом с полосой прокрутки) – линейка ошибок. Если в строке найдена синтаксическая ошибка, то
имеем красный маркер на линейке ошибок; если ошибка некритичная (warning), то цвет маркера желтый. Код может быть предварительно проверен (это чем-то напоминает этап компиляции) - с помощью пункта меню Tools –
Check Code with M-Lint (вызывается специальное окно с перечнем найденных ошибок). Помимо синтаксических есть ошибки времени выполнения. Ошибки времени выполнения выявить более сложно,
33 потому что локальная рабочая область m-функции оказывается потерянной, если ошибка приводит к возврату в рабочую область системы MATLAB. Чтобы определить причину такой ошибки, можно использовать один из следующих приемов:
реализовать вывод результатов промежуточных вычислений на дисплей, удалив в соответствующих операторах точки с запятой, которые подавляют вывод на экран промежуточных результатов;
добавить в m-файл команды keyboard, которые останавливают выполнение m-файла и разрешают проверить и изменить переменные рабочей области вызываемой m-функции. В этом режиме появляется специальное приглашение K». Возврат к выполнению функции реализуется командой return;
закомментировать заголовок функции и выполнить m-файл как сценарий. Это позволяет проследить результаты промежуточных вычислений в рабочей области системы;
использовать отладчик системы МАТLАВ (см. несколько пиктограмм на панели справа) – только в случае отсутствия входных аргументов.
Пример:
1.
Построить и оформить графики трех функций
2 3
3 2
t
1
t
1
)
t
(
z
,
t sin t
1
t ln t
)
t
(
y
, w(t)=z(t)-y(t) на [0;3]. Из серии набранных команд собрать М-файл-сценарий.
Создадим анонимные функции
>> funy=@(x)
(x.^2).*log(x+eps)+(1+x.^2).^(1/3);
>> funz=@(x) (1+x.^3)./(1+x.^2); funw=@(x) funz(x)-funy(x);
Проверим правильность задания
>> X=[0 1 2];[funy(X); funz(X); funw(X)], clear('X') ans =
34 1.0000 1.2599 4.4826 1.0000 1.0000 1.8000 0 -0.2599 -2.6826
Строим графики простейшим путем
>> figure(1);hold on;fplot(funy,[0 3],'- y');fplot(funz,[0 3],'--r');fplot(funw,[0 3],':g');hold off
С помощью контекстного меню (Create…) в History Commands формируем М-файл текущей сессии (лишние команды можно затем стереть в редакторе) (рисунок 10)
Рисунок 10 – Формирование файла текущей сессии
Проведем украшение графика (скажем, заменим желтый цвет линии на синий и добавим легенду) (рисунок 11). Исполним пункт
File/Generate M-file. Откроется новый файл Untitled-2. Скопируем все секции оттуда (кроме первой, начиная со строки %Create Figure).
Сохраним файл под именем MyGraph.m в текущий каталог.
Рисунок 11 – Итоговый график
35
К сожалению, несмотря на полезные конструкции сгенерированного кода, он был исполнен в виде функции, поэтому придется провести ряд изменений, чтобы сделать его работоспособным. В первой секции заменим figure1 = figure() на set(gcf, - ведь полотно графиков уже создано. Во второй секции первая строка приобретет вид set(gca,'Color',[0.502 1 1]);. Также следует убрать третью ось в команде прорисовки осей axis. С секциями «% Create plot» сложнее, здесь придется операцию создания кривой заменять на поиск указанной кривой и установку ее свойств
(см . findobj). Это позволит нам заменить все «пустые» ссылки
(дескрипторы).
Общий вид кода и смасштабированного результата его выполнения (по команде MyGraph в консоли) см. на рисунке 12 снизу и слева.
Рисунок 12 – Общий вид кода
36 2.
Записать в виде М-функции вычисления по формуле
Герона в файле geroncore. Вычислить площадь при наборах (1,2,1.5) и
(1,2,4). function
S=geroncore (a,b,c) p=(a+b+c)/2;
S=sqrt(p.*(p-a).*(p-b).*(p-c));S=real(S); end
>> geroncore ([1 1],[2 2],[1.5 4]) ans =
0.7262 0 3.
Записать в виде М-функции MyTriangle решатель треугольника, т.е. по трем заданным его элементам находящий остальные – число аргументов и результатов переменное.
Приготовим шаблон М-файла, разобьем его на секции (см. иконку «%%+» на второй панели редактора) и сразу его запишем function
[MyS, varargout]=MyTriangle (varagrin)
% Эта функция "решает" треугольник. Стандартный вызов предусматривает
% 4 входных и 4 выходных аргумента. Первые три аргумента - длины трех
% сторон, четвертый - строка вида 'abA', где большими буквами указаны
% углы, а малыми - стороны (угол А противолежит стороне а). Пятый аргумент
% предусматривает дополнительный расчет биссектрис, медиан или высот.
%%
% Основные внутренние переменные a=1;b=1;c=1;A=60;B=A;C=A;MyS=0;MyLegend=
'abc'
;
%%
%Проверка данных
%%
%Запись данных
%%
37
%Вычисление результата
%%
%Вывод результатов varargout(1)=varargin(1); end
Будем наполнять содержание каждой секции. Начнем предварительную проверку данных. При желании формат данных можно выверять тщательнее. Пока же проверяем число аргументов и формат легенды, которая пишется как массив символов в MyLegend.
Для примера мы сохранили закомментированную «точку останова» keyboard. Особенность переменной varargin состоит в том, что она является массивом ячеек (cell array), типом промежуточным между классической записью и массивом. Поэтому для извлечения данных используется функция char. Проще, но менее привычно, было бы записать varargin {4} вместо varargin(4).
%Проверка данных switch nargin case
{0,1,2} warning(
'Слишком мало аргументов'
), flag=false; case
{6,7,8} warning(
'Слишком много аргументов'
), flag=false; otherwise flag=true; if
(nargin=3) MyLegend=char(varargin(4)); end
;
%keyboard;
end
; if
(flag) return
; end
;
MyLegend=MyLegend(1:3);flag=true;NAngles=0; for s=MyLegend switch s case
{
'a'
,
'b'
,
'c'
} case
{
'A'
,
'B'
,
'C'
}
38
NAngles=NAngles+1; otherwise flag=false; end
; if
(flag) break
; end end
; if
(NAngles==3) error(
'По трем углам нельзя решить треугольник'
), return
; end
; if
(flag) error(
'Неправильный формат легенды данных'
), return
; end
;
Для записи данных используем простой цикл и полезную команду eval.
%Запись данных for
(j=1:3) eval(strcat(MyLegend(j),
'='
,num2str(varargin{j}),
'
;'
)); end
;
Для расчета нам нужно разделить секцию на две части, в первой
– описания функций, а во второй – вызывающий их код.
Заключительная команда вызывает функцию geroncore из одноименного М-файла. Мы также оставили точку останова keyboard; поучительно рассмотреть одинаково называющиеся переменные, имеющие разную область действия – см. во время второго останова селектор Workspace-Stack. В стеке функции b=5, а стеке программы b=4 – интерпретатор не путается в действии принципа «локальное имя закрывает глобальное». Теперь программу можно осмысленно запускать, временно раскомментаривая последнюю строчку – например, >> MyTriangle(3,4,5,'abc').
% Сначала список функций, математически все сводится к 5 случаям
%%
% Теорема синусов. Два угла прилежат стороне.
% Находится сторона у первого угла function res=SIN_side (a,alfa,beta)
39 if
(alfa+beta>180) warning(
'Два тупых угла невозможны'
); res=NaN;
return
; end
; res=a.*sin(deg2rad(beta))./sin(deg2rad(alfa+beta))
; end
% Теорема косинусов. Находится сторона возле угла function res=COS_side (b,c,gamma) a1=c.^2-(b.*sin(deg2rad(gamma)))^2; if
(a1<0) warning(
'Такого треугольника не существует'
); res=NaN;
return
; end
; a2=b+sqrt(a1);a1=b-sqrt(a1); if
(a1<0) warning(
'Найдено два таких треугольника!'
); end
; res=a2; end
% Теорема косинусов. Находится сторона против угла function res=COS_opposite (a,b,gamma) res=sqrt(a.^2+b.^2-
2*a.*b.*cos(deg2rad(gamma))); end
% Теорема косинусов. Находится угол (в градусах)
function res=COS_angle (a,b,c) if
((a'Не выполнено неравенство треугольника'
); res=NaN;
return
; end
; keyboard; res=acosd((a.^2+b.^2-c.^2)./(2*a.*b)); end
%%
% Главный блок расчета switch
MyLegend case
'abc'
C=COS_angle(a,b,c); if
(C==NaN) return
; end
;
40
B=COS_angle(a,c,b); if
(B==NaN) return
; end
;
A=COS_angle(c,b,a); if
(A==NaN) return
; end
; case
'abC'
c=COS_opposite(b,a,C);
if
(c==NaN) return
; end
;
B=COS_angle(a,c,b); if
(B==NaN) return
; end
;
A=180-B-C; case
'abA'
c=COS_side(b,a,A);
if
(c==NaN) return
; end
;
B=COS_angle(a,c,b); if
(B==NaN) return
; end
;
C=180-A-B; case
'aBC'
b=SIN_side(a,C,B);
if
(b==NaN) return
; end
; c=COS_opposite(a,b,C);
if
(c==NaN) return
; end
;
A=180-B-C case
'aAB'
C=180-A-B; if
(C<0) return
; end
; b=SIN_side(a,C,B);
if
(b==NaN) return
; end
; c=COS_opposite(a,b,C);
if
(c==NaN) return
; end
; otherwise disp(
'Будет улучшено в след. версиях'
) end
MyS=geroncore(a,b,c);
%[a b c A B C]
Проведем дополнительные вычисления медиан, высот и т.д., в соответствии в 5-м параметром вызова. Обратите внимание, что подфункция может вызывать ранее определенную подфункцию; также –
несмотря на краткость определений функций, каждая должна занимать три строки.
%%
% Дополнительные расчеты по пятому параметру:
% med - медианы, bis - биссектрисы, hhh - высоты
MyS=geroncore(a,b,c);
% Объявление новых функций function res=median(a,b,gamma) res=COS_opposite(a,0.5*b,gamma); end
41 function res=bisectrix(a,b,gamma) res=2*MyS./((a+b).*sind(gamma/2)); end function res=height(c) res=2*MyS./c; end
% Расчет линий for i=1:3, varargout{i}=i; end
; if
(nargin==5) switch varargin{5} case
'med'
varargout{4}=median(b,a,C);varargout{5}=median(a,b
,C); varargout{6}=median(b,c,A); case
'bis'
varargout{4}=bisectrix(c,b,A);varargout{5}=bisectr ix(c,a,B); varargout{6}=bisectrix(a,b,C); case
'hhh'
varargout{4}=height(a);varargout{5}=height(b); varargout{6}=height(c); otherwise
; end end
;
Последняя секция короткая. Список аргументов varargout, как и varargin, - это массив ячеек, поэтому присвоения стоят в {}, а не в привычных (). При выводе графика оси рисуются после построения графика, в противном случае оси будут масштабироваться по данным.
%%
%Вывод результатов varargout{1}=[a A];varargout{2}=[b
B];varargout{3}=[c C];
% Как последний штрих - вывод изображения
42 if
(nargin==5)&&(strcmp(varargin{5},
'qrh'
)) figure(1);plot([1 a.*cosd(C)+1 b+1 1], [1 a.*sind(C)+1 1 1]); axis([0 (a+b+c)/2 0 (a+b+c)/2]); end
Для считывания результата приходится делать вызов программы не слишком красивым, например,
[x,y,z,t,u,v,w]=MyTriangle(3,4,5,’abc’,’grh’)
. Все числа можно было бы расположить компактно в одной переменной- массиве (в качестве результата).
Скачать 1.64 Mb.