Задание:
1.
Построить график функции, заданной в цилиндрических координатах уравнениями:
2
;
0
,
A
cos r
,
sin
Ar
1
z
2
Составить М-файл-сценарий и запустить его несколько раз. На что похожа кривая? Сделать ее «кометой».
Указание: Использовать функции pol2cart, plot3, comet3.
Параметр А>0 свободный.
2.
Оформить функцию myfun1, имеющую два аргумента
(вектор и строку) и один результат – скаляр. Вычисление скаляра зависит от вида строки (предусмотреть три варианта расчета) – например, может рассчитываться среднее арифметическое.
Протестировать М-файл.
3.
Написать программу, вычисляющую средний балл студента. Пользователь должен в консоли вводить имя студента, название предмета и полученный балл. Первым, в качестве базы, вводится название предмета.
Указание: Приветствуется замена циклов встроенными функциями (find,sort), базу данных объявить глобальной переменной
(global) и реализовать как массив ячеек (cell array).
2. Решение некоторых стандартных математических задач
В практике математического моделирования часто приходится сталкиваться с необходимостью решать численно некоторые стандартные задачи. К их числу можно отнести:
Вычислений значений спецфункций
Задача интерполяции – узнать каково значение функции в
43 некоторой точке интервала, если заданы ее значения в нескольких его точках.
Поиск решения системы линейных уравнений, вычисление собственных векторов и пр.
Поиск корней полинома
0
)
x
(
P
n
и вообще, нелинейного уравнения 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 f
2
(x)). Параметр DisplayLevel='iter' позволяет отображать на экране ход итераций процесса и тем самым «почувствовать» успешность оптимизации. fminsearch(fun,x0,options) – ищет безусловный экстремум
44 функций многих переменных.
Для решения более сложных оптимизационных задач см. приложение MATLAB Optimization Toolbox. gradient(F) – дает градиент функции, заданной таблично. dblquad(fun,xmin,xmax,ymin,ymax)
– двукратный интеграл
Римана на прямоугольнике. fft2(X) – быстрое двумерное Фурье-преобразование. odeset – устанавливает параметры для решения обыкновенных дифференциальных уравнений (ОДУ) ode45 – решение ОДУ методом Рунге-Кутты 4-го порядка точности. dde23 – решение ОДУ с запаздыванием.
Пример:
1.
Составить
М-функцию, содержащую несколько независимых вложенных функций. Параметров вызова два – строковая переменная 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;
45 end
%%
function
M5=M5(x,y)
M5=1; end
%%
switch
(S) case
'1'
; case
'2'
; case
'3'
; otherwise disp([S ‘не найдено такого кода’]); end
M=1; end
2.
Вычислить производную в точке x=(x
1
,x
2
) следующей функции W(x): const b
,
a
,
x bx x
)
b
,
x
(
u
),
bx ax
)(
bx
(
)
a x
(
J
)
b
,
x
(
u
)
x
(
W
b
2 1
2 1
2 1
3
/
2
Здесь 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
46
В последней секции исправим (нижняя строка К – коэффициенты, верхняя – координаты точки): 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.
3.
В пределах одного полотна построить два графика: один – контурный для 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])
47 и получим график (рисунок 13). Неровная белая полоса снизу получается по причине Г(+0)=-∞.
Рисунок 13 – Проверочный график
В ядре MATLAB отсутствует функция по вычислению градиента по заданному указателю функции, что вынуждает нас использовать матрицы. Поскольку градиент берется из данных интерполяции, то нужно сначала построить подробную сетку, взять градиент в ее узлах, затем полученные матрицы
вновь интерполировать по нужным нам узлам; в противном случае точность не будет высокой (рисунок 14).
Рисунок 14 - Вычисление градиента по заданному указателю функции
48
%%
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'
);axi s square
;
M3=0; end
В данном случае вызов >> MyMathSolutions('3',[2 3;2 1]).
Неровность снизу контурного графика исчезает; вектора на втором графике показывают направление возрастания функции. Попробуйте сделать шаг вдвое меньше и рассчитать для других параметров.
4.
Вычислить объем n-мерной сферы радиуса 1. Вывести график зависимости V(n)
Математически задача сводится к вычислению n-кратного интеграла, а в отношении программирования более чем уместна рекурсия, т.е. вызов функцией копии самой себя. Однако, MATLAB
49 позволяет вычислять только одно, дву- и троекратные интегралы, а рекурсия в нем не предусмотрена, поэтому воспользуемся методом
Монте-Карло:
1
)
x
,
x
(
,
1
)
x
(
1
)
x
,
x
(
,
0
)
x
(
f
1
)
x
(
)
x
(
плотностью с
величина я
многомерна случайная
)
(
f dx
)
x
(
)
x
(
f
2
)
n
(
V
3 4
)
3
(
V
)
2
(
V
dx dx
)
n
(
V
n k
2
k
R
n
1
x n
1
Заметим, что метод Монте-Карло обладает плохой сходимостью, пропорциональной корню из числа испытаний √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;
%если захотим увидеть ход процесса, каждый сотый шаг
50 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);
5. а) Вычислить все корни уравнения
0 1
x
6
x
3
б) Вычислить все корни трансцендентного уравнения на отрезке
[-5;5]:
A
)
x ln(
)
x cos(
*
x
2
(для некоторых А).
Примечание: излагается решение без М-файла
Рисунок 15 – Отрезок
A
)
x ln(
)
x cos(
*
x
2
>> format long;roots([8 0 -6 1]) ans =
-0.93969262078591 0.76604444311898 0.17364817766693
51
Интересующий нас раздел помощи – 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:
52 the current x satisfies the termination criteria using OPTIONS.TolX of 1.000000e-004 x =
1.46955240774355 fx =
0.08155506454236
Если бы уравнение действительно имело корень, то ордината fx была бы отрицательной. Таким образом, впечатление от графика (о наличии корня) ложное.
6. а) Графически решить систему уравнений (a,b – параметры):
a by x
b y
ax x
2 2
2
б) Найти все решения системы:
1
)
1
z
)(
b y
)(
a x
(
z by ax z
by ax
2 2
2
По пункту (а) – мы используем контурные графики для линий уровня, соответствующих нулю. Добавим в наш М-файл секцию: 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'
);
53 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;
54 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 1
2
z x
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
55 4-Аналитическое решение #(x,y,z): -1.7071 -
2.6909e-006 -1.7071 5-Аналитическое решение #(x,y,z): 1.103 -
2.2059 -3.3089 7.
Решить систему обыкновенных дифференциальных уравнений и начертить траекторию движения частицы (параметр а изменяется от -1 до 1 с шагом 0.5):
t
0 3
2 1
z y
x yt ax dt
/
dz t
sin x
az dt
/
dy t
cos z
ay dt
/
dx
0
t
Особенность решения ОДУ состоит в том, что
выходные вектора многокомпонентны, и точки расположены неравномерно
(последняя проблема решается, однако, командой deval). По параметру а мы создаем семейство решений (рисунок 16) и записываем их в 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,
56
[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),Dat a(:,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)]);
57
Рисунок 16 – Семейство решений
Задание:
1.
Попробуйте п.4. Примера пересчитать с большей точностью, а результаты пересчета отразить на одном графике (для этого немного измените второй фрагмент кода). Поясните полученные результаты. Улучшите код и постройте единый график, адекватный истинным результатам (отобразить их отдельной кривой).
Они определяются аналитическим выражением через гамма- функцию:
1 2
/
n
/
)
n
(
V
2
/
n
2.
Создать М-файл MyMath.m, вычисляющий значение функции
5 3
2
x
5
x
4
x
3
x
2 1
)
x
(
P
, а также дающий его корни.
На отрезке [-1;1] вычислить коэффициенты быстрого Фурье- разложения P(x) (см. fft) кол-ве 10 шт. Сравнить их с коэффициентами классического Фурье-разложения c n
Замечание: Для функции f(x) коэффициенты ряда Фурье задаются формулами (для отрезка [a;b]):
58 2
n
2
n n
b a
n b
a n
b a
5 0
c
,
dx
)
kx a
b
2
sin(
)
x
(
f a
b
2
b
,
dx
)
kx a
b
2
cos(
)
x
(
f a
b
2
a
3.
В рамках того же М-файла реализовать вычисление лапласиана от функции:
))
ay zx
(sin(
zL
z
*
b
)
bz y
x
(
erf
)
xy
(
P
a
)
z
,
y
,
x
(
U
)
xz a
(
z
)
x
1
(
yz x
)
z
,
y
,
x
(
U
)
0
(
5 2
2 2
2 1
Здесь erf – функция ошибок, L
5
– полином Лежандра 5-й степени, a,b- параметры.
Указание: вычислить лапласиан для U
1
аналитически, но не программировать ее.
Сверить результаты аналитики и программирования. Затем адаптировать процедуру для U
2
. Раздел справки (а также математические определения) см. в MATLAB->
Functions -- By Category-> Mathematics->Data Analysis-> Finite
Differences and Integration, а также Mathematics-> Specialized Math.
4.
В рамках того же М-файла реализовать построение 4-х графиков на одном полотне:
1,2,3D-мерного и одного параметрического. На базе данных от функции U
1,2
(x,y,z).
5.
Графическим путем (через пересечение двух кривых) с точностью до 4-го знака найти все корни уравнения 0.1P(x)=U
1
(x,1,1)
– при a=2. На отрезке [-1;2] аппроксимировать U
1
(x,1,1) полиномом 5- й степени Q(x). Найти корни уравнения 0.1P(x)=Q(x).
6.
Пользуясь результатами п.6 Примера, найти все локальные экстремумы функции U
2
(x,y,z) при наборах: (a,b)=(1,1),(1,0),(0,1),(-
1,0).
7. a) Решить следующую систему на отрезке [-5;5] при a=π, b=2,c=-3:
59
2
t
0
,
1 0
1
z y
x
)
cy bx az cos(
dt
/
dz
)
cx bz ay cos(
dt
/
dy
)
cz by ax cos(
dt
/
dx
0
t b) Найти точки Пуанкаре сечения траекторией плоскости z=0.3.
Соединить эти точки с соседними с помощью функции gplot.
Скачать 1.64 Mb.