Файл: М.А. Тынкевич Система Matlab Справочное пособие к курсу Численные методы анализа.pdf
ВУЗ: Не указан
Категория: Не указан
Дисциплина: Не указана
Добавлен: 01.06.2024
Просмотров: 106
Скачиваний: 0
30
8.5. Нули и экстремумы функций
Как известно, решение уравнений и поиск экстремумов функций
– родственные задачи. Так решение системы fi(X)=0, i=1,...,n можно
заменить поиском минимума F( X ) = ∑n fi 2 ( X ) , заведомо равного ну-
i =1
лю, если система имеет решение. С другой стороны, поиск экстремумов функции можно свести к решению системы уравнений относительно нулевых значений ее производных.
Для некоторых классов функций имеются вполне универсальные методы решения: так для полиномов поиск корней реализуется без каких-то условий обобщенным методом Ньютона - функцией roots(...) . В общем случае при великом многообразии методов решение подобных задач отнюдь нетроивиально. В системе реализованы функции лишь для простейших задач без гарантий получения решения.
xmin=fmin(‘имя’,a,b), xmin=fmin(‘имя’,a,b,<опции>), xmin= fmin (‘имя’, a,b,<опции>, p1,...,p10), [xmin,options]=fmin(...) – поиск минимума функции одной переменной на интервале [a,b]; в списке <опции> можно указать option(1)=1 (вывод промежуточных результатов), option(2)- погрешность итераций для аргумента (разница между смежными приближениями; по умолчанию 10-4) и option(14)- максимальное число итераций (по умолчанию 500). Используется метод золотого сечения и параболической интерполяции.
function r=m1(x) |
|
|
|
r=x^2-x-3; |
|
|
|
» xm=fmin(‘m1’,-2,3,[1, 1e-6]) |
|
||
Func |
evals x |
f(x) |
Procedure |
1 |
-0.0901699 |
-2.9017 |
initial |
2 |
1.09017 |
-2.9017 |
golden |
3 |
1.81966 |
-1.5085 |
golden |
4 |
0.5 |
-3.25 |
parabolic |
5 |
0.5 |
-3.25 |
parabolic |
6 |
0.5 |
-3.25 |
parabolic |
xm = |
|
|
|
0.5000
xmin=fmins(‘имя’,x0), xmin=fmins(‘имя’,x0,<опции>), xmin= fmins (‘имя’, x0,<опции> или [], p1,...,p10), [xmin,options]=fmins(...) –
поиск минимума функции нескольких переменных: x0 – начальное приближение; p1,...,p10 – дополнительные параметры; в списке <опции> можно указать option(1)=1 (вывод промежуточных результа-
31
тов), option(2)- погрешность итераций для аргумента (разница между смежными приближениями; по умолчанию 10-4), option(3)- итерационная погрешность для функции и option(14)- максимальное число итераций (по умолчанию 200×n).
function r=m2(x) |
|
|
|
|
|
|
|
|
||
r=(x(1)-1)^2+(x(2)-3)^2; |
|
|
|
|
|
|
|
|||
» [xmin,opt]=fmins(‘m2’,[0.5 2.5]) |
|
|
|
|
|
|
||||
xmin = |
|
|
|
|
|
|
|
|
|
|
1.0000 |
3.0000 |
|
|
|
|
|
|
|
|
|
opt = |
|
|
|
|
|
|
|
|
|
|
0 |
0.0001 |
0.0001 |
0.0000 |
0 |
0 |
0 |
0.0000 |
0 |
|
|
|
72 |
0 |
0 |
0 |
400 |
0 |
0.0000 |
0.1000 |
0 |
|
(возвращаемые здесь опции 8 и 10 определяют минимум функции и число итераций). Используется метод Нелдера-Мида (строится симплекс из n+1 вершины в n –мерном пространстве, берется новая точка внутри или вблизи симплекса и может заменить одну из вершин; процесс повторяется до малого диаметра симплекса).
fzero(‘имя’,x0), fzero(‘имя’,x0,eps), fzero (‘имя’, x0, eps , trace), fzero (‘имя’, x0, eps, trace,p1,...,p10) – поиск действительных корней функции одной переменной при начальном приближении х0 (можно взять и в форме [a,b] при условии f(a)×f(b)<0) с заданной относительной погрешностью; trace=1 – вывод промежуточных результатов. Используется метод дихотомии, хорд и обратной квадратической интерполяции. При поиске корней полинома см. roots.
Процесс поиска корня виден из примера:
» fzero(‘sin’,6,[],1) |
|
|
|
|
|
Func evals |
x |
f(x) |
Procedure |
||
|
1 |
6 |
-0.279415 |
|
initial |
|
2 |
5.83029 |
-0.437567 |
|
search |
|
3 |
6.16971 |
-0.113236 |
|
search |
|
4 |
5.76 |
-0.499642 |
|
search |
|
5 |
6.24 |
-0.0431719 |
search |
|
|
6 |
5.66059 |
-0.583146 |
|
search |
|
7 |
6.33941 |
0.0561963 |
search |
|
Looking for a zero in the interval [5.6606, 6.3394] |
|||||
|
8 |
6.27974 |
-0.00344052 |
interpolation |
|
|
9 |
6.28319 |
1.70244e-006 |
interpolation |
|
|
10 |
6.28319 |
-3.35667e-012 |
interpolation |
|
|
11 |
6.28319 |
-2.44921e-016 interpolation |
||
|
12 |
6.28319 |
2.41961e-015 |
interpolation |
|
ans = |
6.28318530717959 |
|
|
||
» 2*pi |
|
|
|
|
|
ans = |
6.28318530717959 |
|
|
||
32
Заметим, что в случае двух переменных существенную помощь для выбора начального приближения может оказать функция gradient (cм. 8.2) и contour(...), которая будет рассмотрена при ознакомлении с функциями графики.
8.6.Обыкновенные дифференциальные уравнения
Всистеме MatLab предусмотрены специальные средства решения задачи Коши для системобыкновенных дифференциальных урав-
нений, заданных как в явной форме dxdt
= F( t, x ), так и в неявной
M dxdt = F( t,x ) , где М- матрица, - т.н. решатель ОДУ (solver ODE),
обеспечивающий пользователю возможность выбора метода, задания начальных условий и др.
В простейшем варианте достаточно воспользоваться командой
[T,X]=solver(‘F’, [DT], X0, ...), где DT - диапазон интегрирования, X0
– вектор начальных значений, F – имя функции вычисления правых частей системы, solver – имя используемой функции (ode45 - метод Рунге-Кутта 4 и 5-го порядков, ode23 – тот же метод 2 и 3-го порядков, ode113 – метод Адамса - для нежестких систем, ode23s, ode15s – для жестких систем и др.). Версии решателя различаются используемыми методами (по умолчанию относительная погрешность 10-3 и абсолютная 10-6) и соответственно временем и успешностью решения. Под жесткостью здесь понимается повышенное требование к точности – использование минимального шага во всей области интегрирования. При отсутствии информации о жесткости рекомендуется попытаться получить решение посредством ode45 и затем ode15s. Если диапазон DT задан начальным и конечным значением [to , tk], то количество элементов в массиве Т (и в массиве решений Х) определяется необходимым для обеспечения точности шагом; при задании DT в виде [to , t1, t2, ... , tk] или [to :∆t: tk] - указанными значениями.
Например, в простейшем варианте решение уравнения dxdt = t e−t в интервале t [0, 0.5] с начальным условием x(t=0)=1 :
|
function |
f =odu1(t,x) |
|
|
|
|
|
|
|
|
f=t*exp(-t); |
|
|
|
|
|
|
|
» [T,X]=ode45 (‘odu1’, [0, 0.5], 1) |
|
|
|
|
|||
T = |
0 |
0.0125 |
0.0250 |
0.0375 |
0.0500 |
0.0625 |
0.0750 |
0.0875 |
|
0.1000 |
0.1125 |
0.1250 |
0.1375 |
0.1500 |
0.1625 |
0.1750 |
… |
X = |
1.0000 |
1.0000 |
1.0003 |
1.0006 |
1.0012 |
1.0018 |
1.0027 |
1.0036 |
|
1.0046 |
1.0058 |
1.0072 |
1.0086 |
1.0102 |
1.0118 |
1.0136 |
…. |
33
Для иллюстрации решения системы и ряда нестандартных возможностей рассмотрим задачу выравнивания цен по уровню актива в следующей постановке.
Предположим, что изменение уровня актива у пропорционально разности между предложением s и спросом р, т.е. y’=k(s-d), k>0, и что изменение цены z пропорционально отклонению актива у от некоторого уровня y0, т.е. z’=-m(y-y0), m>0. Естественно, что предложение и спрос зависят от цены, например, s(z)=az+s0, d(z)= d0 - cz. Соответственно возникает система дифференциальных уравнений
y’ = k (s(z)-d(z))
z’ = - m (y-y0) .
Вычисление правых частей оформляем функцией:
function f=odu2(t,X) |
|
||
y=X(1); |
z=X(2); |
|
|
a=20; c=10; |
|
s0=10; d0=50; k=0.3; m=0.1; |
|
s=a*z+s0; d= |
d0-c*z; |
y0=19; |
|
f(1)=k*(s-d) |
; f(2)=-m*(y-y0); |
||
f=f'; |
|
% вектор-столбец |
|
Если выполнить решение при y0= 19, z0=2
»[T,Y]=ode45('odu2', [0:0.3:9],[19 2]);
»[T Y]
»plot(T,Y)
будет выведена таблица значений искомых функций: ans = 0 19.0000 2.0000
0.3000 20.7757 1.9731
0.6000 22.4101 1.8949
0.9000 23.7686 1.7714
1.2000 24.7427 1.6126
1.5000 25.2569 1.4313 ...
и их графики (рис. 8.6).
Эта система уравнений относится к числу т.н. автономных (или динамических), ибо независимая переменная в нее явно не входит; соответственно может быть установлена связь между найденными решениями: в параметрическом задании линия y=y(t), z=z(t) определяет фазовую кривую (траекторию) системы - гладкую кривую без самопересечений, замкнутую кривую или точку, которая позволяет судить об устойчивости системы.
Так, установив опции к построению двумерного фазового портрета (функция odephas2)и номера переменных состояния :
»opt=odeset('OutputSel',[1 2], 'OutputFcn','odephas2');
»[T,Y]=ode45('odu2', [0:0.3:9],[19 2],opt);