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

Этот пример показывает, как использовать 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 может эффективно решить уравнение Ван дер Поля. The 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')

Figure contains an axes. The axes with title van der Pol Equation, \mu = 1 contains an object of type line.

Для больших величин μ, задача становится жесткой. Эта метка предназначена для проблем, которые сопротивляются попыткам быть оценены обычными методами. Вместо этого для быстрого интегрирования необходимы специальные числовые методы. The 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')

Figure contains an axes. The axes with title van der Pol Equation, \mu = 1000 contains an object of type line.

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

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 по мере уточнения решения.

The bvpinit функция собирает начальное предположение в форме, которую можно передать решателю bvp4c. Для mesh [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

Figure contains an axes. The axes contains an object of type line.

Это конкретное краевое значение задача имеет ровно два решения. Вы можете получить другое решение, изменив начальные предположения на 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

Figure contains an axes. The axes contains 2 objects of type line. These objects represent Solution 1, Solution 2.

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

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

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

Figure contains an axes. The axes with title An example of Wille and Baker DDE with Constant Delays contains 3 objects of type line. These objects represent y_1, y_2, y_3.

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

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 (при котором требуется моментальный снимок решения). Решите задачу с помощью mesh из 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');

Figure contains an axes. The axes contains an object of type surface.

См. также

| |

Похожие темы