Этот пример показывает, как использовать MATLAB®, чтобы сформулировать и решить несколько различных типов дифференциальных уравнений. MATLAB предлагает несколько числовых алгоритмов, чтобы решить большое разнообразие дифференциальных уравнений:
Задачи с начальными значениями
Краевые задачи
Дифференциальные уравнения с запаздывающим аргументом
Дифференциальные уравнения с частными производными
vanderpoldemo
является функцией, которая задает уравнение Ван дер Поля
type vanderpoldemo
function dydt = vanderpoldemo(t,y,Mu) %VANDERPOLDEMO Defines the van der Pol equation for ODEDEMO. % Copyright 1984-2014 The MathWorks, Inc. dydt = [y(2); Mu*(1-y(1)^2)*y(2)-y(1)];
Уравнение записано как система двух обыкновенных дифференциальных уравнений первого порядка (ОДУ). Эти уравнения оценены для различных значений параметра. Для более быстрого интегрирования необходимо выбрать соответствующий решатель на основе значения.
Поскольку, любой из решателей ОДУ MATLAB может решить уравнение Ван дер Поля эффективно. Решатель ode45
является одним таким примером. Уравнение решено в области с начальными условиями и.
tspan = [0 20]; y0 = [2; 0]; Mu = 1; ode = @(t,y) vanderpoldemo(t,y,Mu); [t,y] = ode45(ode, tspan, y0); % Plot solution plot(t,y(:,1)) xlabel('t') ylabel('solution y') title('van der Pol Equation, \mu = 1')
Для больших значений проблема становится жесткой. Эта метка для проблем, которые сопротивляются попыткам, которые будут оценены с обычными методами. Вместо этого специальные численные методы необходимы для интеграции FAST. ode15s
, ode23s
, ode23t
и функции ode23tb
могут решить жесткие проблемы эффективно.
Это решение уравнения Ван дер Поля для использования ode15s
с теми же начальными условиями. Необходимо протянуть отрезок времени решительно к смочь видеть периодическое перемещение решения.
tspan = [0, 3000]; y0 = [2; 0]; Mu = 1000; ode = @(t,y) vanderpoldemo(t,y,Mu); [t,y] = ode15s(ode, tspan, y0); plot(t,y(:,1)) title('van der Pol Equation, \mu = 1000') axis([0 3000 -3 3]) xlabel('t') ylabel('solution y')
bvp4c
и bvp5c
решают краевые задачи для обыкновенных дифференциальных уравнений.
Функции, взятой в качестве примера, twoode
записали дифференциальное уравнение как систему двух ОДУ первого порядка. Дифференциальное уравнение
type twoode
function dydx = twoode(x,y) %TWOODE Evaluate the differential equations for TWOBVP. % % See also TWOBC, TWOBVP. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. dydx = [ y(2); -abs(y(1)) ];
Функциональный twobc
имеет граничные условия для проблемы: и.
type twobc
function res = twobc(ya,yb) %TWOBC Evaluate the residual in the boundary conditions for TWOBVP. % % See also TWOODE, TWOBVP. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. res = [ ya(1); yb(1) + 2 ];
До вызова bvp4c
необходимо обеспечить предположение для решения, которое вы хотите представленный в mesh. Решатель затем адаптирует mesh, когда это совершенствовало решение.
Функция bvpinit
собирает исходное предположение в форме, можно передать решателю bvp4c
. Для сетки [0 1 2 3 4]
и постоянных предположений и, вызов bvpinit
:
solinit = bvpinit([0 1 2 3 4],[1; 0]);
С этим исходным предположением можно решить проблему с bvp4c
. Оцените решение, возвращенное bvp4c
в некоторых точках с помощью deval
, и затем постройте график получившихся значений.
sol = bvp4c(@twoode, @twobc, solinit); xint = linspace(0, 4, 50); yint = deval(sol, xint); plot(xint, yint(1,:)); xlabel('x') ylabel('solution y') hold on
Эта конкретная краевая задача имеет точно два решения. Можно получить другое решение путем изменения исходных предположений на и.
solinit = bvpinit([0 1 2 3 4],[-1; 0]); sol = bvp4c(@twoode,@twobc,solinit); xint = linspace(0,4,50); yint = deval(sol,xint); plot(xint,yint(1,:)); legend('Solution 1','Solution 2') hold off
dde23
, ddesd
и ddensd
решают дифференциальные уравнения с запаздывающим аргументом с различными задержками. Примеры ddex1
, ddex2
, ddex3
, ddex4
и ddex5
формируют мини-пример при использовании этих решателей.
Пример ddex1
показывает, как решить систему дифференциальных уравнений
Можно представлять эти уравнения с анонимной функцией
ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];
История проблемы (для) является постоянной:
Можно представлять историю как вектор из единиц.
ddex1hist = ones(3,1);
Двухэлементный вектор представляет задержки системы уравнений.
lags = [1 0.2];
Передайте функцию, задержки, историю решения и интервал intergation к решателю как входные параметры. Решатель производит непрерывное решение на целом интервале интегрирования, которое подходит для графического изображения.
sol = dde23(ddex1fun, lags, ddex1hist, [0 5]); plot(sol.x,sol.y); title({'An example of Wille and Baker', 'DDE with Constant Delays'}); xlabel('time t'); ylabel('solution y'); legend('y_1','y_2','y_3','Location','NorthWest');
pdepe
решает дифференциальные уравнения с частными производными в одной пространственной переменной и время. Примеры pdex1
, pdex2
, pdex3
, pdex4
и pdex5
формируют мини-пример при использовании pdepe
.
Эта проблема в качестве примера использует функции pdex1pde
, pdex1ic
и pdex1bc
.
pdex1pde
задает дифференциальное уравнение
type pdex1pde
function [c,f,s] = pdex1pde(x,t,u,DuDx) %PDEX1PDE Evaluate the differential equations components for the PDEX1 problem. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. c = pi^2; f = DuDx; s = 0;
pdex1ic
настраивает начальное условие
type pdex1ic
function u0 = pdex1ic(x) %PDEX1IC Evaluate the initial conditions for the problem coded in PDEX1. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. u0 = sin(pi*x);
pdex1bc
настраивает граничные условия
type pdex1bc
function [pl,ql,pr,qr] = pdex1bc(xl,ul,xr,ur,t) %PDEX1BC Evaluate the boundary conditions for the problem coded in PDEX1. % % See also PDEPE, PDEX1. % Lawrence F. Shampine and Jacek Kierzenka % Copyright 1984-2014 The MathWorks, Inc. pl = ul; ql = 0; pr = pi * exp(-t); qr = 1;
pdepe
требует пространственной дискретизации x
и вектор времен t
(в котором вы хотите снимок состояния решения). Решите проблему с помощью сетки 20 узлов и запросите решение в пяти значениях t
. Извлеките и постройте график первого компонента решения.
x = linspace(0,1,20); t = [0 0.5 1 1.5 2]; sol = pdepe(0,@pdex1pde,@pdex1ic,@pdex1bc,x,t); u1 = sol(:,:,1); surf(x,t,u1); xlabel('x'); ylabel('t'); zlabel('u');