Решение нестандартных ОДУ

Эта страница содержит два примера решения нежестких обыкновенных дифференциальных уравнений с помощью ode45. MATLAB® имеет три решателя для нежестких ОДУ.

  • ode45

  • ode23

  • ode113

Для большинства нежестких проблем ode45 выполняет лучше всего. Однако ode23 рекомендуется для проблем, которые разрешают немного более грубый ошибочный допуск или в присутствии умеренной жесткости. Аналогично, ode113 может быть более эффективным, чем ode45 для проблем со строгими ошибочными допусками.

Если нежесткие решатели занимают много времени, чтобы решить проблему или последовательно приводить интегрирование к сбою, то проблема может быть жесткой. Смотрите Решают Жесткие ОДУ для получения дополнительной информации.

Пример: Нежесткое уравнение Ван дер Поля

Уравнение Ван дер Поля является ОДУ второго порядка

где скалярный параметр. Перепишите это уравнение как систему ОДУ первого порядка путем создания замены. Получившаяся система ОДУ первого порядка

Система ОДУ должна быть закодирована в функциональный файл, который может использовать решатель ОДУ. Общая функциональная подпись функции ОДУ

  dydt = odefun(t,y)

Таким образом, функция должна принять и t и y как входные параметры, даже если это не использует t ни для каких вычислений.

Функциональный файл vdp1.m кодирует использование уравнения Ван дер Поля. Переменные и представлены y(1) и y(2) и двухэлементным вектором - столбцом, dydt содержит выражения для и.

function dydt = vdp1(t,y)
%VDP1  Evaluate the van der Pol ODEs for mu = 1
%
%   See also ODE113, ODE23, ODE45.

%   Jacek Kierzenka and Lawrence F. Shampine
%   Copyright 1984-2014 The MathWorks, Inc.

dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];

Решите ОДУ с помощью функции ode45 на временном интервале [0 20] с начальными значениями [2 0]. Вывод является вектором - столбцом моментов времени t и массив решения y. Каждая строка в y соответствует времени, возвращенному в соответствующей строке t. Первый столбец y соответствует, и второй столбец к.

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

Постройте график решений для и против t.

plot(t,y(:,1),'-o',t,y(:,2),'-o')
title('Solution of van der Pol Equation (\mu = 1) using ODE45');
xlabel('Time t');
ylabel('Solution y');
legend('y_1','y_2')

Функция vdpode решает ту же проблему, но это принимает пользовательское заданное значение для. Уравнения Ван дер Поля становятся жесткими как увеличения. Например, со значением необходимо использовать жесткий решатель, такой как ode15s, чтобы решить систему.

Пример: Нежесткие эйлеровы уравнения

Эйлеровы уравнения для твердого тела без внешних сил являются стандартной тестовой задачей для решателей ОДУ, предназначенных для нежестких проблем.

Уравнения

Функциональный файл rigidode задает и решает эту систему уравнений первого порядка по временному интервалу [0 12], с помощью вектора начальных условий [0; 1; 1], соответствующий начальным значениям, и. Локальная функция f(t,y) кодирует систему уравнений.

rigidode вызывает ode45 без выходных аргументов, таким образом, решатель использует выходную функцию по умолчанию odeplot, чтобы автоматически построить график точек решения после каждого шага.

function rigidode
%RIGIDODE  Euler equations of a rigid body without external forces.
%   A standard test problem for non-stiff solvers proposed by Krogh.  The
%   analytical solutions are Jacobian elliptic functions, accessible in
%   MATLAB.  The interval here is about 1.5 periods; it is that for which
%   solutions are plotted on p. 243 of Shampine and Gordon.
%
%   L. F. Shampine and M. K. Gordon, Computer Solution of Ordinary
%   Differential Equations, W.H. Freeman & Co., 1975.
%
%   See also ODE45, ODE23, ODE113, FUNCTION_HANDLE.

%   Mark W. Reichelt and Lawrence F. Shampine, 3-23-94, 4-19-94
%   Copyright 1984-2014 The MathWorks, Inc.

tspan = [0 12];
y0 = [0; 1; 1];

% solve the problem using ODE45
figure;
ode45(@f,tspan,y0);

% --------------------------------------------------------------------------

function dydt = f(t,y)
dydt = [    y(2)*y(3)
   -y(1)*y(3)
   -0.51*y(1)*y(2) ];


Решите нежесткие Эйлеровы уравнения путем вызывания функции rigidode.

rigidode
title('Solution of Rigid Body w/o External Forces using ODE45')
legend('y_1','y_2','y_3','Location','Best')

Смотрите также

| |

Похожие темы

Была ли эта тема полезной?