Лабораторное занятие №11- №15
Проведение вычислений с помощью М-файлов
1. Создание и отладка М-файлов.
MATLAB является интерпретирующим языком непосредственных вычислений, т.е. выражения, которые вы вводите, интерпретируются и вычисляются. Поскольку все вычисления в MATLAB выполняются с двойной точностью, формат вывода может управляться с помощью следующих команд
После вызова одного из приведенных выше форматов он сохраняется до вызова другого. Команда format compact подавляет большинство пустых строк, позволяя большее количество информации вывести на экран или страницу. Она не зависит от других команд формата. Команда eval чрезвычайна полезна, позволяя обмениваться между процедурами (особенно, внешними) строкой-параметром, которая принуждает систему к выполнению заданной этой строкой команды (или серии команд). Напомним, что MATLAB различает большие и маленькие буквы в именах команд, функций и переменных, а спецзнаки могут быть записаны через обратный слэш.
Для простых операций удобен интерактивный режим, но если вычисления нужно многократно выполнять или необходимо реализовывать сложные алгоритмы, то следует использовать m-файлы MATLAB (расширение файла состоит из одной буквы m). Познакомимся со script-m-файлами (или сценариями)– текстовыми файлами, содержащими инструкции на языке MATLAB, подлежащими исполнению в автоматическом пакетном режиме. Создать такой файл удобнее с помощью редактора системы MATLAB (не забудем способ создания, исходя из окна истории команд). Он вызывается из командного окна системы MATLAB командой меню File/New/M-file (или самой левой кнопкой на полосе инструментов, на которой изображен чистый белый лист бумаги). Записанные в script-файлы команды будут выполнены, если в командной строке ввести имя script-файла (без расширения). Переменные, определяемые в командном окне и переменные, определяемые в сценариях, составляют единое рабочее пространство системы MATLAB, причем переменные, определяемые в сценариях, являются глобальными, их значения заместят значения таких же переменных, которые были использованы до вызова данного script-файла. Текст реально имеющихся m-файлов (системы MATLAB, например polar, или ваших собственных) можно просмотреть с помощью команды type <имя_функции>.
Однако, большинство m-файлов являются файлами-функциями, т.е. программами. При том же порядке вызова в отличие от простого линейного скрипта, хотя бы тот и содержал, например, циклические конструкции, нормальный m-файл содержит тело, окаймленное следующими декларациями:
function [out1, out2, ...] = funname(in1, in2, ...)
end
Тело может быть представлено в таком формате:
=
%Строки основного комментария, доступного по команде help
% Далее три четыре основных повторяющихся типа конструкций
% Объявление собственных переменных и команды MATLAB
<переменная=выражение>; <команда>
% Объявление подфункций (пишется обычно с отступом)
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, 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 (вызывается специальное окно с перечнем найденных ошибок). Помимо синтаксических есть ошибки времени выполнения. Ошибки времени выполнения выявить более сложно, потому что локальная рабочая область m-функции оказывается потерянной, если ошибка приводит к возврату в рабочую область системы MATLAB. Чтобы определить причину такой ошибки, можно использовать один из следующих приемов:
реализовать вывод результатов промежуточных вычислений на дисплей, удалив в соответствующих операторах точки с запятой, которые подавляют вывод на экран промежуточных результатов;
добавить в m-файл команды keyboard, которые останавливают выполнение m-файла и разрешают проверить и изменить переменные рабочей области вызываемой m-функции. В этом режиме появляется специальное приглашение K». Возврат к выполнению функции реализуется командой return;
закомментировать заголовок функции и выполнить m-файл как сценарий. Это позволяет проследить результаты промежуточных вычислений в рабочей области системы;
использовать отладчик системы МАТLАВ (см. несколько пиктограмм на панели справа) – только в случае отсутствия входных аргументов.
Пример:
Построить и оформить графики трех функций  , 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 =
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 формируем М-файл текущей сессии (лишние команды можно затем стереть в редакторе)
Проведем украшение графика (скажем, заменим желтый цвет линии на синий и добавим легенду). Исполним пункт File/Generate M-file. Откроется новый файл Untitled-2. Скопируем все секции оттуда (кроме первой, начиная со строки %Create Figure). Сохраним файл под именем MyGraph.m в текущий каталог.
К сожалению, несмотря на полезные конструкции сгенерированного кода, он был исполнен в виде функции, поэтому придется провести ряд изменений, чтобы сделать его работоспособным. В первой секции заменим figure1 = figure() на set(gcf, - ведь полотно графиков уже создано. Во второй секции первая строка приобретет вид set(gca,'Color',[0.502 1 1]);. Также следует убрать третью ось в команде прорисовки осей axis. С секциями «% Create plot» сложнее, здесь придется операцию создания кривой заменять на поиск указанной кривой и установку ее свойств (см . findobj). Это позволит нам заменить все «пустые» ссылки (дескрипторы).
Общий вид кода и смасштабированного результата его выполнения (по команде MyGraph в консоли) см. на рис. снизу и слева.
Записать в виде М-функции вычисления по формуле Герона в файле 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
Записать в виде М-функции 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';
%%
%Проверка данных
%%
%Запись данных
%%
%Вычисление результата
%%
%Вывод результатов
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'}
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)
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
warning('Не выполнено неравенство треугольника');
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;
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
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}=bisectrix(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];
% Как последний штрих - вывод изображения
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’). Все числа можно было бы расположить компактно в одной переменной-массиве (в качестве результата).
Задание:
Построить график функции, заданной в цилиндрических координатах уравнениями:  . Составить М-файл-сценарий и запустить его несколько раз. На что похожа кривая? Сделать ее «кометой».
Указание: Использовать функции pol2cart, plot3, comet3. Параметр А>0 свободный.
Оформить функцию myfun1, имеющую два аргумента (вектор и строку) и один результат – скаляр. Вычисление скаляра зависит от вида строки (предусмотреть три варианта расчета) – например, может рассчитываться среднее арифметическое. Протестировать М-файл.
Написать программу, вычисляющую средний балл студента. Пользователь должен в консоли вводить имя студента, название предмета и полученный балл. Первым, в качестве базы, вводится название предмета.
Указание: Приветствуется замена циклов встроенными функциями (find,sort), базу данных объявить глобальной переменной (global) и реализовать как массив ячеек (cell array).
2. Решение некоторых стандартных математических задач
В практике математического моделирования часто приходится сталкиваться с необходимостью решать численно некоторые стандартные задачи. К их числу можно отнести:
Вычислений значений спецфункций
Задача интерполяции – узнать каково значение функции в некоторой точке интервала, если заданы ее значения в нескольких его точках.
Поиск решения системы линейных уравнений, вычисление собственных векторов и пр.
Поиск корней полинома Pn(x)=0 и вообще, нелинейного уравнения f(x)=0
Вычисление коэффициентов разложения Фурье (в т.ч. быстрое Фурье–преобразование)
Вычисление значений производных, дивергенций, определенных интегралов от функций
Задачи поиска наименьшего и наибольшего значений (т.е. локальных и глобальных экстремумов)
Поиск решений обыкновенных дифференциальных уравнений и их систем (задача Коши)
Краевые задачи для некоторых уравнений в частных производных
Всю необходимую информацию можно почерпнуть в разделе справки MATLAB-Mathematics, а также в разделе MATLAB–Functions_by_Category–Mathematics. Кратко о некоторых полезных функциях (по их имени легко восстановить формат и особенности вызова):
legendre(n,X) – вычисление значений полиномов Лежандра порядка n в точках вектора Х.
yy = spline(x,Y,xx) – интерполяция кубическими сплайнами yy(xx) (с возможной экстраполяцией), хх – как скаляр, так и вектор.
poly(A) – дает коэффициенты характеристического многочлена матрицы А.
roots(p) – дает корни полинома, заданного вектором коэффициентов.
fzero(fun,x0) – находит нуль функции, заданной указателем fun, в окрестности x0.
optimset – устанавливает параметры оптимизации (например, для fzero – заметим, что решение алгебраических уравнений может быть сведено к оптимизационной задаче; например, f(x*)=0↔min f2(x)). Параметр DisplayLevel='iter' позволяет отображать на экране ход итераций процесса и тем самым «почувствовать» успешность оптимизации.
fminsearch(fun,x0,options) – ищет безусловный экстремум функций многих переменных.
Для решения более сложных оптимизационных задач см. приложение MATLAB Optimization Toolbox.
gradient(F) – дает градиент функции, заданной таблично.
dblquad(fun,xmin,xmax,ymin,ymax) – двукратный интеграл Римана на прямоугольнике.
fft2(X) – быстрое двумерное Фурье-преобразование.
odeset – устанавливает параметры для решения обыкновенных дифференциальных уравнений (ОДУ)
ode45 – решение ОДУ методом Рунге-Кутты 4-го порядка точности.
dde23 – решение ОДУ с запаздыванием.
Пример:
Составить М-функцию, содержащую несколько независимых вложенных функций. Параметров вызова два – строковая переменная S, по которой определяется вызываемая подфункция, и массив коэффициентов К.
function M=MyMathSolutions(S,K)
%%
function M1=M1(x,y)
M1=1;
end
%%
function M2=M2(x,y)
M2=2;
end
%%
function M3=M3(x,y)
M3=3;
end
%%
function M4=M4(x,y)
M4=4;
end
%%
function M5=M5(x,y)
M5=1;
end
%%
switch (S)
case '1'
;
case '2'
;
case '3'
;
otherwise
disp([S ‘не найдено такого кода’]);
end
M=1;
end
Вычислить производную в точке x=(x1,x2) следующей функции W(x):
Здесь J – функция Бесселя первого рода, Г – гамма-функция, u – определенная пользователем функция
Прочитаем справку по используемым спецфункциям. Выпишем первую секцию (введем «защиту от дураков» в виде модуля, поскольку Г(x<0) не определено, как и в нуле):
%%
function M1=M1(x,p)
user_def_fun=@(x,b) x(1)+b*x(2)+norm(x).^b;
M1=user_def_fun(x,p(2)).*besselj(2/3,x(1)-p(1))-gamma(abs(p(2)*x(2)))*(p*x');
end
В последней секции исправим (нижняя строка К – коэффициенты, верхняя – координаты точки):
case '1'
M=M1(K(1,1:2),K(2,1:2));
....%M=1;
После вызова-теста >> MyMathSolutions('1',[1 1;1 –1]) получим верный ответ 0. Пробный расчет дает MyMathSolutions('1',[2 3;2 1])= –14.
В пределах одного полотна построить два графика: один – контурный для W(x), второй – график скоростей для градиента gradW(x). Принять a=1.5,b=2.3, площадка – первая четверть круга радиуса π в центром в (0,0). Шаг сетки для графика скоростей равен π/4.
Для простоты можно было бы использовать циклы, поскольку фактическая трехмерность матрицы данных мешает использовать ezcontourf. И все-таки имеет место такая короткая формулировка (данные транспонируются):
%%
function M2=M2(x,y)
if (x.^2+y.^2>pi^2) M2=NaN;
else M2=real(M1([x.' y.'],K(2,1:2)));
end
end
И в первом приближении нужно добавить в последнюю секцию:
case '2'
ezcontourf(@M2,[0,pi,0,pi]);M=2;
В качестве проверки запустим >> MyMathSolutions('2',[2 3;2 1]) и получим график. Неровная белая полоса снизу получается по причине Г(+0)=-∞.
В ядре MATLAB отсутствует функция по вычислению градиента по заданному указателю функции, ч то вынуждает нас использовать матрицы. Поскольку градиент берется из данных интерполяции, то нужно сначала построить подробную сетку, взять градиент в ее узлах, затем полученные матрицы вновь интерполировать по нужным нам узлам; в противном случае точность не будет высокой.
%%
function M3=M3(x,y)
 [X,Y]=meshgrid(0:0.1:pi);Z=X;
for i=1:length(X),
for j=1:length(Y),
Z(i,j)=real(M2(X(i,j),Y(i,j)));
end
end
[Zx,Zy]=gradient(Z,0.05,0.05);
xbase=0:pi/4:pi;[xx,yy]=meshgrid(xbase);
Fx = griddata(X,Y,Zx,xx,yy);Fy = griddata(X,Y,Zy,xx,yy);
scrsz = [1 1 0.4 0.8].*get(0,'ScreenSize');figure('Name',...
'Simulation Plot Window','NumberTitle',...
'off','Position',scrsz);
subplot(2,1,1);contourf(X,Y,Z);colorbar;axis square;
subplot(2,1,2);quiver(xx,yy,Fx,Fy,'Color','r');axis square;
M3=0;
end
В данном случае вызов >> MyMathSolutions('3',[2 3;2 1]). Неровность снизу контурного графика исчезает; вектора на втором графике показывают направление возрастания функции. Попробуйте сделать шаг вдвое меньше и рассчитать для других параметров.
Вычислить объем n-мерной сферы радиуса 1. Вывести график зависимости V(n).
Математически задача сводится к вычислению n-кратного интеграла, а в отношении программирования более чем уместна рекурсия, т.е. вызов функцией копии самой себя. Однако, MATLAB позволяет вычислять только одно, дву- и троекратные интегралы, а рекурсия в нем не предусмотрена, поэтому воспользуемся методом Монте-Карло:
Заметим, что метод Монте-Карло обладает плохой сходимостью, пропорциональной корню из числа испытаний √N.
Соответствующие фрагменты выглядят так (была проведена нормировка на число π):
function M4=MySphereVolume(x,p,flag)
% x- вещественный аргумент, p - целые параметры, например, размерность
% при начальном вызове x равен радиусу, а p - скаляр
MassOfPoint=@(x) 1;
function Region=Region(x)
if (sum(x.*x)>1) Region=0; else Region=MassOfPoint(x); end
end
N=2;summa=1;medsumma=0;medsumma1=0.5;eps=1e-8;
while ((abs(medsumma-medsumma1)>eps)&&(N*sqrt(eps)<1))
ksi=rand(1,p);summa=summa+Region(ksi);N=N+1;
medsumma=medsumma1;medsumma1=summa/N;
%если захотим увидеть ход процесса, каждый сотый шаг
if (flag&&(~mod(N,100))) ['процессинг: ' num2str(medsumma) ...
' ' num2str(medsumma1)], end;
end
%вторым аргументом идет погрешность расчета
M4=[medsumma1*((2*x).^p) 1/sqrt(N)];
end
….
case '4'
Z=[];
for n=1:10, res=MySphereVolume(1,n,false), Z=[Z res(1)/pi]; end;
plot(Z);
а) Вычислить все корни уравнения 8x3-6x+1=0.
б) Вычислить все корни трансцендентного уравнения на отрезке [-5;5]: x*cos(x)+ln(x2)=A (для некоторых А).
Примечание: излагается решение без М-файла
>> format long;roots([8 0 -6 1])
ans =
-0.93969262078591
0.76604444311898
0.17364817766693
Интересующий нас раздел помощи – MATLAB->Mathematics->Function Functions->Minimizing Functions and Finding Zeros. Чтобы «почувствовать» задачу «на лету» построим график:
A=1; transcen=@(x) x.*cos(x)+log(x.^2)-A; x=-5:0.01:5;plot(x,transcen(x)), grid on
Можно попробовать вручную, мышкой, найти корень. Для этого в графическом редакторе выбрать пиктограмму  , в контекстном меню на кривой выбрать Selection Style – Mouse Position, Display Style – DataTip, выбрать стартовую точку, и, не отпуская кнопку мыши, провести по кривой. График позволяет нам выделить участки перемены знака. Найдем третий по возрастанию корень:
>> fzero(transcen,[4 5])
ans =
4.25053115278136
Однако график не позволяет нам определить, является ли второй корень двойным, или мы имеем два близких корня. Здесь, несмотря на многочисленность параметров (см. optimset) fzero, нам придется использовать функцию оптимизации fminbnd:
>> transcen1=@(x) -transcen(x);[x,fx]=fminbnd(transcen1,1,2,optimset('Display','iter'))
Func-count x f(x) Procedure
1 1.38197 0.0935767 initial
2 1.61803 0.11398 golden
3 1.23607 0.170065 golden
4 1.47297 0.0815729 parabolic
5 1.47129 0.0815597 parabolic
6 1.46961 0.0815551 parabolic
7 1.46955 0.0815551 parabolic
8 1.46952 0.0815551 parabolic
Optimization terminated:
the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004
x =
1.46955240774355
fx =
0.08155506454236
Если бы уравнение действительно имело корень, то ордината fx была бы отрицательной. Таким образом, впечатление от графика (о наличии корня) ложное.
а) Графически решить систему уравнений (a,b – параметры):
б) Найти все решения системы:
По пункту (а) – мы используем контурные графики для линий уровня, соответствующих нулю. Добавим в наш М-файл секцию:
function M5=M5(a,b)
%Решение графически системы двух уравнений
%Сетка адаптирована к параметрам a,b
f1=@(x,y) x.^2+y.^2-a*x-b^2;f2=@(x,y) x+b*y-1;
x=fix(-b+a/2-1):b/100:fix(b+a/2+1);y=fix(-b-1):0.001:fix(b+1);
[X,Y]=meshgrid(x,y);Z1=f1(X,Y);Z2=f2(X,Y);
figure(1);hold on;contour(X,Y,Z1,[0 0],'-r','DisplayName','1');
contour(X,Y,Z2,[0 0],'-g','DisplayName','2');
axis equal; grid on; legend show; hold off;
M5=0;
end
…
case '5'
M5(K(1,1),K(2,1));
По пункту (б) – задача сводится к поиску нулевого экстремума. Добавим еще секцию:
%%
function M6=M6(a,b)
%Решение системы некольких уравнений
%Область поиска локальных экстремумов разбивается на
%несколько случайным числом
f1=@(x) a*x(1)+ b*x(2)-x(3);
f2=@(x) a*x(1).^2+b*x(2).^2-x(3).^2;
f3=@(x) (x(1)+a).*(x(2)+b).*(x(3)+1)-1;
%Критерий корня - экстремум суммы квадратов
OurFun=@(x) sqrt(f1(x).^2+f2(x).^2+f3(x).^2);
Res=[];
for n=1:100,
display(['Номер итерации: ' int2str(n)]), pause(0.05);
opt=optimset('Display','off');p=10*rand(1,3)-5;
[x,fx]=fminsearch(OurFun,p,opt);
if (abs(fx)<1e-4) Res=[Res;x fx]; end;
end;
%убрать повторы
if (size(Res,1)==0) display('Решений нет'),
else
Res=Res';Rest=Res(:,1);%keyboard;
for p=Res,
flag=true;
for q=Rest,
%keyboard;
if (norm(p-q)<1e-3) flag=false; end;
end;
if (flag) Rest=[Rest p]; end;
end;
end;
%проверка правых частей, 3 последних строки
Rest=[Rest(1:3,:);Rest(1:3,:)];n=0;
for p=Rest, p(4)=f1(p(1:3));p(5)=f2(p(1:3));p(6)=f3(p(1:3));
n=n+1;Rest(:,n)=p;
display([int2str(n) '-Аналитическое решение #(x,y,z): ' num2str(p(1:3)')]),
%display(['Погрешности от уравнений: ' num2str(p(4:6)')]),
end;
M6=Rest;
End
…
case '6'
M=M6(K(1,1),K(2,1));
При а=1 и b=2 вызов W=MyMathSolutions('6',K) дает всего 5 решений (решения №3,4 можно получить аналитически – при y=0 x=z=±2-1/2-1):
…
Номер итерации: 100
1-Аналитическое решение #(x,y,z): -0.93181 1.8636 2.7954
2-Аналитическое решение #(x,y,z): 0.16215 -0.32434 -0.4865
3-Аналитическое решение #(x,y,z): -0.29292 1.8644e-005 -0.29288
4-Аналитическое решение #(x,y,z): -1.7071 -2.6909e-006 -1.7071
5-Аналитическое решение #(x,y,z): 1.103 -2.2059 -3.3089
Решить систему обыкновенных дифференциальных уравнений и начертить траекторию движения частицы (параметр а изменяется от -1 до 1 с шагом 0.5):
Особенность решения ОДУ состоит в том, что выходные вектора многокомпонентны, и точки расположены неравномерно (последняя проблема решается, однако, командой deval). По параметру а мы создаем семейство решений и записываем их в 3-х-мерную матрицу. Входная матрица К равна [0 1;3 2].
%%
function M7=M7(a,b)
function RateX=RateX(x,y,z,t)
RateX=a*y-z.*cos(t);
end
function RateY=RateY(x,y,z,t)
RateY=a*z-x.*sin(t);
end
function RateZ=RateZ(x,y,z,t)
RateZ=a*x-y.*t;
end
function Rate=Rate(time,arg)
Rate=[RateX(arg(1),arg(2),arg(3),time);...
RateY(arg(1),arg(2),arg(3),time); RateZ(arg(1),arg(2),arg(3),time)];
end
Data=[];Ti=linspace(0,pi,100)';
for a=-1:0.5:1,
[T,Y]=ode45(@Rate,[0 pi],b);
%Шаг интегрирования нельзя задавать постоянным.
%Интерполируем данные
Y1 = interp1(T,Y(:,1),Ti,'spline'); Y2 = interp1(T,Y(:,2),Ti,'spline');
Y3 = interp1(T,Y(:,3),Ti,'spline'); Data=cat(3,Data,[Y1 Y2 Y3]);
end;
%Данные подсчитаны. Приступим к построению графиков на одном полотне
scrsz = get(0,'ScreenSize');
figure('Name','Траектории','Position',[0.1*scrsz(3:4) 0.8*scrsz(3:4)]);
for j=1:5,
subplot(2,3,j);
plot(Ti,Data(:,1,j),'-',Ti,Data(:,2,j),'-.', Ti,Data(:,3,j),':');
grid on;legend('X(t)','Y(t)','Z(t)');legend('show');
title(gca,['Параметр равен: a=' num2str(j/2-1.5)]);
end;
% На лишнем месте построим 3D-кривую для а=0.5
subplot(2,3,6);h=plot3(Data(:,1,4),Data(:,2,4),Data(:,3,4));
set(h,'Color','m'); axis square;grid on;box('on');title('Параметр 0.5');
legend({'Траектория'},'Orientation','horizontal','location','NorthEast');
M7=1;
end
…
case '7'
M=M7(K(1,1),[K(1,2);K(2,2);K(2,1)]);
Задание:
Попробуйте п.4. Примера пересчитать с большей точностью, а результаты пересчета отразить на одном графике (для этого немного измените второй фрагмент кода). Поясните полученные результаты. Улучшите код и постройте единый график, адекватный истинным результатам (отобразить их отдельной кривой). Они определяются аналитическим выражением через гамма-функцию:  .
Создать М-файл MyMath.m, вычисляющий значение функции P(x)=1+2x+3x2+4x3+5x5, а также дающий его корни. На отрезке [-1;1] вычислить коэффициенты быстрого Фурье-разложения P(x) (см. fft) кол-ве 10 шт. Сравнить их с коэффициентами классического Фурье-разложения cn.
Замечание: Для функции f(x) коэффициенты ряда Фурье задаются формулами (для отрезка [a;b]):
В рамках того же М-файла реализовать вычисление лапласиана от функции:
Здесь erf – функция ошибок, L5 – полином Лежандра 5-й степени, a,b- параметры.
Указание: вычислить лапласиан для U1 аналитически, но не программировать ее. Сверить результаты аналитики и программирования. Затем адаптировать процедуру для U2. Раздел справки (а также математические определения) см. в MATLAB-> Functions -- By Category-> Mathematics->Data Analysis-> Finite Differences and Integration, а также Mathematics-> Specialized Math.
В рамках того же М-файла реализовать построение 4-х графиков на одном полотне: 1,2,3D-мерного и одного параметрического. На базе данных от функции U1,2(x,y,z).
Графическим путем (через пересечение двух кривых) с точностью до 4-го знака найти все корни уравнения 0.1P(x)=U1(x,1,1) – при a=2. На отрезке [-1;2] аппроксимировать U1(x,1,1) полиномом 5-й степени Q(x). Найти корни уравнения 0.1P(x)=Q(x).
Пользуясь результатами п.6 Примера, найти все локальные экстремумы функции U2(x,y,z) при наборах: (a,b)=(1,1),(1,0),(0,1),(-1,0).
a) Решить следующую систему на отрезке [-5;5] при a=π, b=2,c=-3:
b) Найти точки Пуанкаре сечения траекторией плоскости z=0.3. Соединить эти точки с соседними с помощью функции gplot.
|
