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

Этот пример показывает, как использовать 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. ode15s, ode23s, ode23t и функции 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

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

dde23, ddesd и 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');

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

| |

Похожие темы