Этот пример показывает, как использовать 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 может эффективно решить уравнение Ван дер Поля. The 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')
Для больших величин , задача становится жесткой. Эта метка предназначена для проблем, которые сопротивляются попыткам быть оценены обычными методами. Вместо этого для быстрого интегрирования необходимы специальные числовые методы. The 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 по мере уточнения решения.
The bvpinit
функция собирает начальное предположение в форме, которую можно передать решателю bvp4c
. Для mesh [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
сформировать мини- руководство на использовании этих решателей.
The ddex1
пример показывает, как решить систему дифференциальных уравнений
Можно представить эти уравнения с помощью анонимной функции
ddex1fun = @(t,y,Z) [Z(1,1); Z(1,1)+Z(2,2); y(2)];
История задачи (для ) является постоянным:
Можно представлять историю как вектор таковых.
ddex1hist = ones(3,1);
Двухэлементный вектор представляет задержки в системе уравнений.
lags = [1 0.2];
Передайте функцию, задержки, историю решений и интервал интегрирования в решатель как входы. Решатель вырабатывает непрерывное решение на протяжении всего интервала интегрирования, которое подходит для графического изображения.
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
(при котором требуется моментальный снимок решения). Решите задачу с помощью 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');