Дифференциальные уравнения

В этом примере показано, как использовать MATLAB®, чтобы сформулировать и решить несколько различных типов дифференциальных уравнений. MATLAB предлагает несколько числовых алгоритмов, чтобы решить большое разнообразие дифференциальных уравнений:

  • Задачи с начальными значениями

  • Краевые задачи

  • Дифференциальные уравнения с запаздывающим аргументом

  • Дифференциальные уравнения с частными производными

Задача с начальными значениями

vanderpoldemo функция, которая задает уравнение Ван дер Поля

d2ydt2-μ(1-y2)dydt+y=0.

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)];

Уравнение записано как система двух обыкновенных дифференциальных уравнений первого порядка (ОДУ). Эти уравнения оценены для различных значений параметра μ. Для более быстрого интегрирования необходимо выбрать соответствующий решатель на основе значения μ.

Для μ=1, любой из решателей ОДУ MATLAB может решить уравнение Ван дер Поля эффективно. ode45 решатель является одним таким примером. Уравнение решено в области [0,20] с начальными условиями y(0)=2 и dydt|t=0=0.

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. ode15sode23sode23t, и ode23tb функции могут решить жесткие задачи эффективно.

Это решение уравнения Ван дер Поля для μ=1000 использование ode15s с теми же начальными условиями. Необходимо протянуть отрезок времени решительно к [0,3000] смочь видеть периодическое перемещение решения.

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 записали дифференциальное уравнение как систему двух ОДУ первого порядка. Дифференциальное уравнение

d2ydt2+|y|=0.

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 имеет граничные условия для проблемы: y(0)=0 и y(4)=-2.

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] и постоянные предположения y(x)=1 и y'(x)=0, вызов 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

Эта конкретная краевая задача имеет точно два решения. Можно получить другое решение путем изменения исходных предположений в y(x)=-1 и y'(x)=0.

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

Дифференциальные уравнения с запаздывающим аргументом

dde23ddesd, и ddensd решите дифференциальные уравнения с запаздывающим аргументом с различными задержками. Примеры ddex1, ddex2, ddex3, ddex4, и ddex5 сформируйте мини-пример при использовании этих решателей.

ddex1 пример показывает, как решить систему дифференциальных уравнений

y1'(t)=y1(t-1)y2'(t)=y1(t-1)+y2(t-0.2)y3'(t)=y2(t).

Можно представлять эти уравнения анонимной функцией

ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];

История проблемы (для t0) является постоянным:

y1(t)=1y2(t)=1y3(t)=1.

Можно представлять историю как вектор из единиц.

ddex1hist = ones(3,1);

Двухэлементный вектор представляет задержки системы уравнений.

lags = [1 0.2];

Передайте функцию, задержки, историю решения и интервал интегрирования [0,5] к решателю как входные параметры. Решатель производит непрерывное решение на целом интервале интегрирования, которое подходит для графического вывода.

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 определяет дифференциальное уравнение

π2ut=x(ux).

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 настраивает начальное условие

u(x,0)=sinπx.

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 настраивает граничные условия

u(0,t)=0,

πe-t+xu(1,t)=0.

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');

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

| |

Похожие темы