Как решать дифференциальные уравнения в матлабе

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

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

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

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

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

d 2 y d t 2 — μ ( 1 — y 2 ) d y d t + 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 и dy dt | 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 записали дифференциальное уравнение как систему двух ОДУ первого порядка. Дифференциальное уравнение

d 2 y d t 2 + | 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 пример показывает, как решить систему дифференциальных уравнений

y 1 ‘ ( t ) = y 1 ( t — 1 ) y 2 ‘ ( t ) = y 1 ( t — 1 ) + y 2 ( t — 0 . 2 ) y 3 ‘ ( t ) = y 2 ( t ) .

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

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

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

y 1 ( t ) = 1 y 2 ( t ) = 1 y 3 ( 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 определяет дифференциальное уравнение

π 2 ∂ u ∂ t = ∂ ∂ x ( ∂ u ∂ x ) .

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

π e — t + ∂ ∂ x u ( 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');

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

Похожие темы

• Выбор решателя ОДУ
• Решение Краевых задач
• Решение дифференциальных уравнений с частными производными

Открытый пример

У вас есть модифицированная версия этого примера. Вы хотите открыть этот пример со своими редактированиями?

Документация MATLAB

Поддержка

• MATLAB Answers
• Помощь в установке
• Отчеты об ошибках
• Требования к продукту
• Загрузка программного обеспечения

© 1994-2021 The MathWorks, Inc.

• Условия использования
• Патенты
• Торговые марки
• Список благодарностей

Для просмотра документации необходимо авторизоваться на сайте
Войти
Памятка переводчика

1. Если смысл перевода понятен, то лучше оставьте как есть и не придирайтесь к словам, синонимам и тому подобному. О вкусах не спорим.

2. Не дополняйте перевод комментариями “от себя”. В исправлении не должно появляться дополнительных смыслов и комментариев, отсутствующих в оригинале. Такие правки не получится интегрировать в алгоритме автоматического перевода.

3. Сохраняйте структуру оригинального текста — например, не разбивайте одно предложение на два.

4. Не имеет смысла однотипное исправление перевода какого-то термина во всех предложениях. Исправляйте только в одном месте. Когда Вашу правку одобрят, это исправление будет алгоритмически распространено и на другие части документации.

5. По иным вопросам, например если надо исправить заблокированное для перевода слово, обратитесь к редакторам через форму технической поддержки.

Документация

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

В partial differential equation (УЧП) функция, решаемая для, зависит от нескольких переменных, и дифференциальное уравнение может включать частные производные, взятые относительно каждой из переменных. Дифференциальные уравнения с частными производными полезны для моделирования волн, теплового потока, жидкой дисперсии и других явлений с пространственным поведением, которое изменяется в зависимости от времени.

Какие типы УЧП можно решить с MATLAB ?

MATLAB ® Решатель УЧП pdepe решает начальные краевые задачи для систем УЧП в одной пространственной переменной x и время t . Можно думать о них как об ОДУ одной переменной, которые также изменяются относительно времени.

pdepe использует неофициальную классификацию для 1D уравнений, которые она решает:

• Уравнениями с производной времени является parabolic. Примером является уравнение тепла ∂ u ∂ t = ∂ 2 u ∂ x 2 .
• Уравнениями без производной времени является elliptic. Примером является уравнение Лапласа ∂ 2 u ∂ x 2 = 0 .

pdepe требует по крайней мере одного параболического уравнения в системе. Другими словами, по крайней мере одно уравнение в системе должно включать производную времени.

pdepe также решает определенные 2D и 3-D задачи, которые уменьшают до 1D проблем из-за угловой симметрии (см. описание аргумента для симметрии постоянный m для получения дополнительной информации.

Partial Differential Equation Toolbox™ расширяет эту функциональность к обобщенным проблемам в 2D и 3-D с Дирихле и Неймановыми граничными условиями.

Решение 1D УЧП

1D УЧП включает функциональный u (x, t) , который зависит от времени t и одна пространственная переменная x. Решатель УЧП MATLAB pdepe решает системы 1D параболических и эллиптических УЧП формы

c ( x , t , u , ∂ u ∂ x ) ∂ u ∂ t = x − m ∂ ∂ x ( x m f ( x , t , u , ∂ u ∂ x ) ) + s ( x , t , u , ∂ u ∂ x ) .

Уравнение имеет свойства:

• УЧП содержат для t₀ ≤ t ≤ t_f и a ≤ x ≤ b .
• Пространственный интервал [a, b] должен быть конечным.
• m может быть 0, 1, или 2, соответствуя плите , цилиндрической , или сферической симметрии, соответственно. Если m> 0 , то a ≥ 0 должен также содержать.
• Коэффициент f ( x , t , u , ∂ u ∂ x ) термин потока и s ( x , t , u , ∂ u ∂ x ) характеристики выброса.
• Термин потока должен зависеть от частной производной ∂u / ∂ x .

Связь частных производных относительно времени ограничивается умножением диагональной матрицей c ( x , t , u , ∂ u ∂ x ) . Диагональными элементами этой матрицы является или нуль или положительный. Элемент, который является нулем, соответствует эллиптическому уравнению, и любой другой элемент соответствует параболическому уравнению. Должно быть по крайней мере одно параболическое уравнение. Элемент c, который соответствует параболическому уравнению, может исчезнуть в изолированных значениях x, если они — точки mesh (точки, где решение оценено). Разрывы в c и s из-за материальных интерфейсов разрешены при условии, что точка mesh помещается в каждый интерфейс.

Процесс решения

Решить УЧП с pdepe , необходимо задать коэффициенты уравнения для c, f, и s, начальных условий, поведения решения на контурах и сетки точек, чтобы оценить решение на. Вызов функции sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan) использование эта информация, чтобы вычислить решение на заданную mesh:

• m постоянная симметрия.
• pdefun определяет решаемые уравнения.
• icfun задает начальные условия.
• bcfun задает граничные условия.
• xmesh вектор из пространственных значений для x.
• tspan вектор из временных стоимостей для t.

Вместе, xmesh и tspan векторы формируют 2D сетку это pdepe оценивает решение на.

Уравнения

Необходимо описать УЧП в стандартной форме, ожидаемой pdepe . Написанный в этой форме, можно прочитать значения коэффициентов c, f и s.

В MATLAB можно закодировать уравнения с функцией формы

function [c,f,s] = pdefun(x,t,u,dudx) c = 1; f = dudx; s = 0; end

В этом случае pdefun определяет уравнение ∂ u ∂ t = ∂ 2 u ∂ x 2 . Если существует несколько уравнений, то c, f и s являются векторами с каждым элементом, соответствующим одному уравнению.

Начальные условия

В начальное время t = t ₀, для всего x, компоненты решения удовлетворяют начальным условиям формы

u ( x , t 0 ) = u 0 ( x ) .

В MATLAB можно закодировать начальные условия с функцией формы

function u0 = icfun(x) u0 = 1; end

В этом случае u0 = 1 задает начальное условие u ₀ (x, t ₀) = 1 . Если существует несколько уравнений, то u0 вектор с каждым элементом, задающим начальное условие одного уравнения.

Граничные условия

За пределами x = a или x = b, для всего t, компоненты решения удовлетворяют граничным условиям формы

p ( x , t , u ) + q ( x , t ) f ( x , t , u , ∂ u ∂ x ) = 0.

q (x, t) является диагональной матрицей с элементами, которые являются или нулем или никогда не обнуляют. Обратите внимание на то, что граничные условия описываются в терминах потока f, а не частная производная u относительно x. Кроме того, этих двух коэффициентов p (x, t, u) и q (x, t) , только p может зависеть от u.

В MATLAB можно закодировать граничные условия с функцией формы

function [pL,qL,pR,qR] = bcfun(xL,uL,xR,uR,t) pL = uL; qL = 0; pR = uR - 1; qR = 0; end 

pL и qL коэффициенты для левого контура, в то время как pR и qR коэффициенты для правильного контура. В этом случае bcfun задает граничные условия

u L ( x L , t ) = 0 u R ( x R , t ) = 1

Если существует несколько уравнений, то выходные параметры pL , qL , pR , и qR векторы с каждым элементом, задающим граничное условие одного уравнения.

Опции интегрирования

Свойства интегрирования по умолчанию в решателе УЧП MATLAB выбраны, чтобы решить типичные проблемы. В некоторых случаях можно улучшать производительность решателя путем переопределения этих значений по умолчанию. Для этого использовать odeset создать options структура. Затем передайте структуру pdepe как последний входной параметр:

sol = pdepe(m,pdefun,icfun,bcfun,xmesh,tspan,options)

Из опций для базового решателя ОДУ ode15s , только показанные в следующей таблице доступны для pdepe .

Оценка решения

После того, как вы решаете уравнение с pdepe , MATLAB возвращает решение как трехмерный массив sol , где sol(i,j,k) содержит k компонент th решения оценен в t(i) и x(j) . В общем случае можно извлечь k компонент решения th с командой u = sol(. k) .

Mesh времени, которую вы задаете, используется просто для выходных целей и не влияет на внутренние временные шаги, взятые решателем. Однако пространственная mesh, которую вы задаете, может влиять на качество и скорость решения. После решения уравнения можно использовать pdeval оценивать структуру решения, возвращенную pdepe с различной пространственной mesh.

Пример: уравнение тепла

Примером параболического УЧП является уравнение тепла в одной размерности:

∂ u ∂ t = ∂ 2 u ∂ x 2 .

Это уравнение описывает рассеяние тепла для 0 ≤ x ≤ L и t ≥ 0 . Цель состоит в том, чтобы решить для температуры u ( x , t ) . Температура является первоначально ненулевой константой, таким образом, начальное условие

Кроме того, температура является нулем на левом контуре, и ненулевой на правильном контуре, таким образом, граничные условия

Чтобы решить это уравнение в MATLAB, необходимо закодировать уравнение, начальные условия и граничные условия, затем выбрать подходящую mesh решения прежде, чем вызвать решатель pdepe . Или можно включать необходимые функции как локальные функции в конце файла, (как в этом примере) или сохранить их как отдельные, именованные файлы в директории на пути MATLAB.

Уравнение кода

Прежде чем можно будет закодировать уравнение, необходимо убедиться что именно в форме pdepe решатель ожидает:

c ( x , t , u , ∂ u ∂ x ) ∂ u ∂ t = x — m ∂ ∂ x ( x m f ( x , t , u , ∂ u ∂ x ) ) + s ( x , t , u , ∂ u ∂ x ) .

В этой форме уравнение тепла

1 ⋅ ∂ u ∂ t = x 0 ∂ ∂ x ( x 0 ∂ u ∂ x ) + 0 .

Таким образом, значения коэффициентов следующие:

• m = 0
• c = 1
• f = ∂ u ∂ x
• s = 0

Значение m передается в качестве аргумента pdepe , в то время как другие коэффициенты закодированы в функции для уравнения, которое является

function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end 

(Примечание: Все функции включены как локальные функции в конце примера.)

Условие начальной буквы кода

Начальная функция условия для уравнения тепла присваивает постоянное значение для u 0 . Эта функция должна принять вход для x , даже если это не использовано.

function u0 = heatic(x) u0 = 0.5; end 

Граничные условия кода

Стандартная форма для граничных условий ожидается pdepe решатель

p ( x , t , u ) + q ( x , t ) f ( x , t , u , ∂ u ∂ x ) = 0 .

Написанный в этой форме, граничные условия для этой проблемы

u ( 0 , t ) + ( 0 ⋅ f ) = 0 ,

( u ( L , t ) — 1 ) + ( 0 ⋅ f ) = 0 .

Так значения для p и q

• p L = u L , q L = 0 .
• p R = u R — 1 , q R = 0 .

Соответствующая функция затем

function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end 

Выберите Solution Mesh

Используйте пространственную сетку 20 точек и сетку времени 30 точек. Поскольку решение быстро достигает устойчивого состояния, моменты времени рядом t = 0 более близко расположены вместе, чтобы получить это поведение в выходе.

L = 1; x = linspace(0,L,20); t = [linspace(0,0.05,20), linspace(0.5,5,10)];

Решите уравнение

Наконец, решите уравнение с помощью симметрии m , уравнение PDE, начальное условие, граничные условия и сетки для x и t .

m = 0; sol = pdepe(m,@heatpde,@heatic,@heatbc,x,t);

Постройте решение

Используйте imagesc визуализировать матрицу решения.

colormap hot imagesc(x,t,sol) colorbar xlabel('Distance x','interpreter','latex') ylabel('Time t','interpreter','latex') title('Heat Equation for $0 \le x \le 1$ and $0 \le t \le 5$','interpreter','latex')

Локальные функции

function [c,f,s] = heatpde(x,t,u,dudx) c = 1; f = dudx; s = 0; end function u0 = heatic(x) u0 = 0.5; end function [pl,ql,pr,qr] = heatbc(xl,ul,xr,ur,t) pl = ul; ql = 0; pr = ur - 1; qr = 0; end

Примеры УЧП и файлы

Несколько доступных файлов в качестве примера служат превосходными начальными точками для наиболее распространенных 1D проблем УЧП. Чтобы исследовать и запустить примеры, используйте приложение Дифференциальных уравнений В качестве примера. Чтобы запустить это приложение, ввести

odeexamples

Чтобы открыть отдельный файл для редактирования, ввести

edit exampleFileName.m

Чтобы запустить пример, ввести

exampleFileName

Эта таблица содержит список доступных файлов УЧП в качестве примера.

Файл в качестве примера

Простой УЧП, который иллюстрирует формулировку, расчет и графический вывод решения.

Проблема, которая включает разрывы.

Проблема, которая требует вычисления значений частной производной.

Система двух УЧП, решение которых имеет пограничные слои и в концах интервала и изменяется быстро для маленького t.

Система УЧП со ступенчатыми функциями как начальные условия.

Ссылки

[1] Skeel, R. D. и М. Берзиньш, «Метод для Пространственной Дискретизации Параболических уравнений в Одной Пространственной переменной», SIAM Journal на Научном и Статистическом Вычислении, Издании 11, 1990, стр 1–32.

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

Похожие темы

• Решите один УЧП
• Решите систему УЧП

Открытый пример

У вас есть модифицированная версия этого примера. Вы хотите открыть этот пример со своими редактированиями?

Решение системы дифференциальных уравнений методом эйлера матлаб

Доброго времени суток! Сегодня мы поговорим о решении ОДУ (обыкновенных дифференциальных уравнений) в Matlab. Перед тем как мы начнём обсуждать данную тему, советую вам ознакомиться с темой: Численное дифференцирование в Matlab, чтобы лучше понимать теоретическую составляющую решения ОДУ.

Обыкновенные дифференциальные уравнения

С помощью дифференциальных уравнений можно описать разные задачи: движения системы, взаимодействующих материальных точек, химической кинетики и т.д. Различают три типа задач для систем диф. уравнений:

• Задача Коши
• Краевая задача
• Задача на собственные значения

Кратко расскажу о их сути:

Задача Коши предполагает дополнительные условия в виде значения функции в определённой точке.
Краевая задача подразумевает поиск решения на заданном отрезке с краевыми (граничными) условиями в концах интервала или на границе области.
Задача на собственные значения — помимо искомых функций и их производных, в уравнение входят дополнительное несколько неизвестных параметров, которые являются собственными значениями.

Видео: Метод Эйлера Скачать

Методы решения дифференциальных уравнений

Решение ОДУ в Matlab и не только, в первую очередь, сводится к выбору порядка численного метода решения. Порядок численного метода не связан с порядком дифференциального уравнения. Высокий порядок у численного метода означает его скорость сходимости.

В случае большого интервала, с помощью алгоритмов с низким порядком сжимают интервал с решениями и находят приблизительные корни, а затем уже уточняют корни с помощью методов с высоким порядком.

Решение обыкновенных дифференциальных уравнений в Matlab можно реализовать «своими ручками», прописав алгоритм по разным схемам. Но также в Matlab есть встроенные функции, выполняющие все стандартные задачи.

Метод Рунге-Кутта первого порядка

Методы Рунге-Кутта представляют собой разложения в ряд Тейлора и от количества использованных элементов ряда зависит порядок этого метода. Следовательно, помимо Рунге-Кутта первого порядка, вы сможете увидеть методы других порядков. Иногда их называют другими именами.

Например, Метод Рунге-Кутта первого порядка, также известен как Метод Эйлера или Метод ломаных. Информацию о его математическом и графическом представлении советую поискать в гугл. Мы же поговорим о том, как Метод Рунге-Кутта первого порядка реализуется в Matlab для решения ОДУ. Например:

Решить и привести график ошибки уравнения y’ = y*x методом Рунге-Кутта первого порядка (Методом Эйлера, Методом ломаных).

Погрешность Метода Рунге-Кутта 1 порядка

» data-medium-file=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?fit=300%2C236&ssl=1″ data-large-file=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?fit=622%2C489&ssl=1″ loading=»lazy» src=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/%D0%A0%D1%83%D0%BD%D0%B3%D0%B5-1-%D0%BF%D0%BE%D0%B3%D1%80%D0%B5%D1%88%D0%BD%D0%BE%D1%81%D1%82%D1%8C.png?resize=622%2C489″ alt=»Погрешность метода 1 порядка» width=»622″ height=»489″ srcset=»https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?w=629&ssl=1 629w, https://i2.wp.com/codetown.ru/wp-content/uploads/2017/02/Рунге-1-погрешность.png?resize=300%2C236&ssl=1 300w» sizes=»(max-width: 622px) 100vw, 622px» data-recalc-dims=»1″ />
На данном графике показана зависимость величины ошибки от шага.

Метод Рунге-Кутта второго порядка

Также известен как Метод Эйлера-Коши. Как видите, во второй части уравнения происходит обращения к следующему шагу. Но как тогда быть, если нам ещё не известен следующий шаг? Всё просто. Метод Рунге-Кутта второго порядка — это всё тот же метод первого порядка, однако, на половине шага происходит нахождение «первичного» решения, а затем происходит его уточнение. Это позволяет поднять порядок скорости сходимости до двух.

Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта второго порядка.

По сравнению с Рунге-Куттом первого порядка изначальная ошибка уже гораздо меньше.

Мы не будем говорить о третьем порядке, потому что задачи на третий порядок встречаются редко, но если будет необходимо, пишите в комментариях, выложу.

Метод Рунге-Кутта четвёртого порядка

Метод Рунге-Кутта четвёртого порядка считается самым распространённым. Тем не менее, работает он аналогично второму и третьему порядку.

Решить и привести график ошибки уравнения u’ = u*x методом Рунге-Кутта четвёртого порядка.

Как видите, на последней картинке размерность ошибки на столько мала, что пришлось воспользоваться loglog() для лучшей видимости.

Решение ОДУ в Matlab стандартными средствами

Стоит отметить, что мы с вами разобрали только один самый известный метод решения ОДУ с разными порядками. Однако, методов очень много.

Для решения дифференциальных уравнений и систем в MATLAB предусмотрены следующие функции:

ode45 (f, interval, X0, [options])
ode23 (f, interval, X0, [options])
ode113 (f, interval, X0, [options])
ode15s (f, interval, X0, [options])
ode23s (f, interval, X0, [options])
ode23t (f, interval, X0, [options])
ode23tb (f, interval, X0, [options])

Входными параметрами этих функций являются:

• f — вектор-функция для вычисления правой части уравнения системы уравнений;
• interval — массив из двух чисел, определяющий интервал интегрирования дифференциального уравнения или системы;
• Х0 — вектор начальных условий системы дифференциальных уравнений;
• options — параметры управления ходом решения дифференциального уравнения или системы.

• массив Т — координаты узлов сетки, в которых ищется решение;
• матрицу X, i-й столбец которой является значением вектор-функции решения в узле Тi.

В функции ode45 реализован метод Рунге-Кутта 4-5 порядка точности, в функции ode23 также реализован метод Рунге-Кутта, но 2-3 порядка, а функция ode113 реализует метод Адамса.

Для решения жёстких систем предназначены функция ode15s, в которой реализован метод Гира, и функция ode23s, реализующая метод Розенброка. Для получения более точного решения жёсткой системы лучше использовать функцию ode15s. Для решения системы с небольшим числом жёсткости можно использовать функцию ode23t, а для грубой оценки подобных систем служит функция ode23tb.

Символьное решение обыкновенных дифференциальных уравнений произвольного порядка осуществляет функция dsolve r = dsolve(‘eq1,eq2,…’, ‘cond1,cond2,…‘, ‘v’)
Пример использования:

На этом мы закончим. Если остались вопросы, задавайте их в комментариях. Также вы можете скачать исходники чтобы лучше понять тему: «Решение ОДУ в Matlab».

Видео: Численное решение системы дифференциальных уравнений(задачи Коши) Скачать

Библиотека функций MATLAB’а

Важным средством, позволяющим ускорить разработку программ, являются библиотеки стандартных программ. Есть такая библиотека и в MATLAB’е.

Надо научиться находить в HELP’е описания нужных программ. В частности в конце этого занятия вам будет предложено написать программу решения дифференциальных уравнений по методу Эйлера. Написав свою программу, вам надо будет найти функцию MATLAB’а, решающую эту же задачу и используя ее сравнить решения полученные вашей программой и библиотечной

Что еще почитать

Как вы понимаете, на одном занятии невозможно освоить принципы отладки и тестирования в полной мере. Если вы намерены совершенствовать свое мастерство программиста, вам придется ознакомиться с более серьезными работами по этому вопросу, чем эти краткие заметки.

Надо сказать, что такой литературы, которая учила бы принципам программирования, а не синтаксическим правилам какого-либо языка программирования не так уж много. «Программистов часто учат программировать и редко – отладке» (ван Тассел).

Могу порекомендовать следующие книги, в которых вопросы отладки и тестирования подробно рассмотрены:

1. Ван Тассел Д. Стиль, разработка, эффективность, отладка и испытания программ. 1981
2. Керниган, Пайк. Практика программирования. 2004
3. Иванова Г.С. Технология программирования. 2002

Задания

Написать программу решения дифференциальных уравнений по методу Эйлера. При отладке и тестировании вести учет ошибок. Получить решение этой же задачи, используя библиотечную функцию MATLAB’а. Сравнить полученные результаты

№	Задача Коши для уравнения первого порядка
Задача Коши для уравнения высшего порядка

Справочная информация

Обыкновенное дифференциальное уравнение 1-го порядка

связывает независимую переменную t,искомую функцию y и ее производную. Решение дифференциальногоуравнениязаключается в отысканиифункции y = y(t), обращающей это уравнение в тождество на конечном или бесконечном интервале (a, b). Различают общее и частные решения дифференциального уравнения. Общее решение имеет вид y = y(t, C), где C — произвольная постоянная интегрирования. Его графическим отображением является семейство кривых (см. рис.1), называемых интегральными. Каждая интегральная кривая является отображением частного решения, соответствующего своему значению постоянной C. Для выделения частного решения из множества общего решения необходимо задать начальное условие

Такая постановка задачи решения дифференциальных уравнений называется задачей Коши.

Для решения задачи Коши существует множество методов, которые делятся на одношаговые и многошаговые. Все они позволяют получить искомое решение в виде таблично заданной функции, в той или иной мере согласующееся с истинным частным решением (см. рис.2). Эти группы методов различаются объемом информации, которая используется для вычисления координат очередной точки табличной функции. Одношаговые методы используют значения функции и ее производной только в одной предыдущей точке, в то время как многошаговые — в нескольких. К одношаговым методам решения задачи Коши относятся метод Эйлера, модифицированный метод Эйлера, методы Рунге-Кутта и другие.

Задача Коши для нормальной системы обыкновенных дифференциальных уравнений записывается в виде

или в матричной форме следующим образом

Система дифференциальных уравнений связывает независимую переменную x, искомые функции y₁, y₂, . y_n и их первые производные. В данном случае решение задачи Коши заключается в отыскании функций y₁ = y₁(x), y₂ = y₂(x). y_n = y_n(x), обращающих каждое уравнение системы в тождество на конечном или бесконечном интервале (a, b) и удовлетворяющих начальным условиям.

Такая форма записи задачи Коши является канонической для систем обыкновенных дифференциальных уравнений. К ней могут быть приведены как любые другие формы представления систем дифференциальных уравнений, разрешенных относительно старших производных, так и аналогичные дифференциальные уравнения высших порядков. Приведение дифференциальных уравнений высших порядков к нормальной системе дифференциальных уравнений осуществляется по следующей схеме. Если дана задача Коши следующего вида

что является задачей Коши для нормальной системы дифференциальных уравнений.

Для решения такой задачи Коши используются те же методы, что для обыкновенных дифференциальных уравнений 1-го порядка. Это обуславливается тем, что матричная форма записи задачи Коши для нормальной системы полностью совпадает с ее формулировкой для этих уравнений. Единственным отличием является то, что вместо функций y(t) и f(t,y) используются вектор-функции y и f, состоящие из n функций y₁(t), y₂(t). y_n(t) и f₁(t,y₁. y_n), f₂(t,y₁. y_n). f_n(t,y₁. y_n), соответственно. При этом расчетные схемы методов и оценки их погрешностей сохраняются.

Он является старейшим методом решения задачи Коши и заключается в последовательном применении следующих формул

геометрическая интерпретация которых представлена на рис.3. В точке t₀ вычисляется значение производной dy/dt через f(t,y), которое определяет угол наклона касательной к графику точного решения задачи Коши. Следующая точка численного решения определяется как точка на этой касательной с абсциссой t₁ = t₀ + h. В компактном виде эти соотношения записываются следующим образом

Метод Эйлера — самый простой метод численного интегрирования. Он относится к методам первого порядка точности, поскольку его решение совпадает с истинным только в том случае, когда истинное решение является линейной функцией y = a₁ + a₂t. Его погрешность e_k на каждом шаге пропорциональна величине h 2 . Это обусловлено тем, что в качестве направления, определяющего положение следующей точки численного решения, используется касательная в крайней левой точке каждого отрезка [t_k,t_k+1]. Из рис.3 видно, что для получения более точного численного решения надо использовать некоторое промежуточное направление между направлениями касательных в крайних точках рассматриваемого отрезка.

Соотношения метода Эйлера для нормальной системы в матричной форме имеют вид

или в развернутой форме

Геометрическая интерпретация работы метода Эйлера решения задачи Коши для нормальной системы идентична его геометрической интерпретации для дифференциальных уравнений 1-го порядка. Однако в данном случае движение осуществляется вдоль некоторойгиперкривойв (n+1)-мерном пространстве переменных t, y₁, y₂. y_n.

Программное обеспечение

Ниже приведен текст подпрограммы на алгоритмическом языке FORTRAN, реализующей метод Эйлера:

Subroutine Eu(n, dt, k, Fun, t, Y)

В качестве параметров в подпрограмме используются:

n — порядок системы дифференциальных уравнений;

dt — длина отрезка интегрирования уравнений;

k — количество шагов на длине отрезка интегрирования;

Fun — имя внешней подпрограммы типа Subroutine, с помощью которой вычисляются значения вектора правой части f(t, y)системы;

t — значение аргумента системы в начале отрезка интегрирования при обращении к подпрограммам и в конце этого отрезка (t + dt)после отработки подпрограмм;

Y — массив (n элементов) значений вектора решения y(t)системы уравнений в начале отрезка интегрирования при обращении к подпрограммам и в конце этого отрезка после отработки подпрограмм.

Подпрограмма Fun оформляется в виде

Subroutine Fun(n, t, Y, F)

Dimension Y(n), F(n)

Здесь F — массив значений вектора правой части f(t, y) системы уравнений.

Subroutine Eu(n, dt, k, Fun, t, Y)

Call Fun(n, t, Y, F1)

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

· [t,Y] = ode23(fun, tspan, y0) – решает задачу Коши для системы обыкновенных дифференциальных уравнений с использованием комбинации методов Рунге-Кутта относительно невысокого порядка (2-го и 3-го);

· [t,Y] = ode45(fun, tspan, y0) — решает задачу Коши для системы обыкновенных дифференциальных уравнений с использованием комбинации методов Рунге-Кутта более высокого порядка (4-го и 5-го).

Входные аргументы указанных функций:

fun — функция для вычисления вектора правой части системы, интерфейс которой должен иметь вид:

где t – независимый параметр системы, Y- вектор значений неизвестных функций, а выходной параметр F – вектор вычисленных компонент правой части;

tspan — вектор, определяющий интервал интегрирования системы. Если этот вектор состоит из двух компонент, то их значения задают начало и конец интервала. Если требуется, чтобы в качестве решения были выданы значения искомых функций в конкретно заданных точках интервала, то число компонент должно быть не менее трех (можно задавать в виде имени массива или при помощи конструктора массивов, например в виде [2:0.1:5]);

y0 — вектор начальных значений.

t — вектор, состоящий из значений независимой переменной, которым соответствуют вычисленные значения решения, помещенные в массив Y (если вектор tspan состоит лишь из двух компонент, то количество значений автоматически определяется функцией с целью наиболее точного отображения на графике);

Y — двумерный массив, каждый столбец которого представляет последовательность вычисленных значений одной из искомых функций.

Параметр fun, являющийся фактическим параметром, указывающим предварительно написанную функцию для вычисления правой части системы, задается в виде:

@fname , где fname — имя m-функции.

Точность получаемого решения по умолчанию регулируется условием обеспечения относительной погрешности не более 10 -3 или абсолютной погрешности не более 10 -6 . Если требуется задать иные требования точности, то функции следует вызывать с дополнительным аргументом options:

[t,Y] = ode23(fun, tspan, y0, options) или

[t,Y] = ode45(fun, tspan, y0, options) .

Значение аргумента options, представляющего собой структурную переменную, содержащую значения управляющих параметров, определяющих режимы работы этих функций, формируется предварительно до вызова функции интегрирования системы при помощи оператора вида

где namei — название управляющего параметра, а valuei — новое присваиваемое ему значение. Среди управляющих параметров имеются RelTol и AbsTol , задающие относительную и абсолютную погрешности решения систем. Так, например, если предполагается, что вызываемая функция должна обеспечить относительную погрешность решения не более 10 -5 или абсолютную не более 10 -8 , то следует выполнить оператор

а затем использовать переменную options в качестве фактического аргумента при решении системы дифференциальных уравнений.

Видео: 5 Численное решение дифференциальных уравнений Part 1 Скачать

Решение систем обыкновенных дифференциальных уравнений в среде MATLAB. Часть 1

В среде MATLAB можно решать системы диффуров с начальными условиями, краевые задачи, а также решать дифференциальные уравнения в частных производных с помощью инструмента PDE toolbox.

В данном обзоре речь пойдет лишь о системах дифференциальных уравнений с начальными условиями, то есть о задаче Коши. В англоязычной литературе это называется Initial Value Problem.

• каким образом записывать системы диффуров
• как задать начальные условия
• временной интервал
• какой получать результат решения для дальнего использования

Решать системы обыкновенных дифференциальных уравнений можно как в MATLAB, так и в Simulink.

В первую очередь, следует определиться, использовать для решения Matlab и его текстовый редактор, или Simulink, где те же системы дифференциальных уравнений могут быть записаны в виде функциональных блоков.

Выбор ваш должен зависеть от задачи. Если Вы, например, хотите смоделировать какой-либо объект управления, описываемый системой диффуров, то в данном случае имеет смысл использовать именно Simulink, так как Вам, впоследствии, понадобиться синтез, например, системы управления, и Simulink подойдет здесь как нельзя лучше.

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

Рассмотрим синтаксис решателей matlab.В качестве аргументов следует подать правую часть системы в виде MATLAB-функции.

На рисунке показан требуемый вид системы, когда выражены старшие производные.

Системы, чей вид отличается от требуемого, следует преобразовать к таковому.

Если функция простая, то её можно записать прямо в поле аргумента, однако, когда речь идёт о системах уравнений, имеет смысл записывать систему уравнений в виде отдельной функции, в том числе и в виде отдельного м-файла. Об этом мы поговорим чуть позже и на конкретном примере.

Также подается интервал времени, на котором будет найдено решение. Интервал задаётся строкой из двух чисел: начальной величины независимого аргумента t и его конечного значения.

Далее задаются начальные условия. Значения всех неизвестных искомых переменных в начале расчёта задаются в виде столбца соответствующей размерности.

Далее, при необходимости, задаются опции. Вот тут и раскрываются широкие возможности MATLAB по настройке решателя. Помимо управления точностью и величиной шага, имеется возможность обрабатывать данные в процессе вычисления, а также выполнять скрипты по завершению вычисления. Но ещё более полезным является опция отслеживания событий по условию, более подробно поговорим об этом дальше. Также есть другие специальные опции, которые могут быть использованы при решении определённых типов систем.

Вы могли заметить, что название функции — odeXY – это обозначение для всех решателей, которых всего 8 штук. В данном ролике мы пользоваться решателем ode45, соответствующего численному по методу Дормана-Принса 4(5). Этого решателя достаточно для подавляющего большинства задач. Остальные решатели будут подробно рассмотрены в приложении к задачам соответствующих типов позже.

Перейдем к примерам.

Рассмотрим 2 примера:

• решение дифференциального уравнения первого порядка.
• решение системы двух дифференциальных уравнений второго порядка.

В качестве уравнение первого порядка рассмотрим логистическое уравнение Ферхюльста, которое описывает динамику численности популяции. Суть уравнения такова: скорость прироста населения y пропорциональна количеству населения, однако лимитирована максимальной численностью популяции.

Забавный факт: Ферхюльст назвал это уравнение логистическим, и никто до сих пор не знает почему, ибо сам Ферхюльст об этом никому не рассказал.

Решение этого дифференциального уравнения выглядит следующим образом.

Пишем функцию в явном виде, задаём интервал расчёта и записываем начальное условие. Пару слов о записи функции подобным образом. Знак собаки в матлабе является оператором создания функции соответствующих переменных. Вы задаёте аргументы функции и саму функцию через пробел, как показано на рисунке.

Перейдем в окно MATLABа и посмотрим, как это выглядит.

Так выглядит скрипт:

Так выглядит график решения дифференциального уравнения:

В качестве примера решения системы, состоящей из двух дифференциальных уравнений второго порядка, рассмотрим шарик, подвешенный на пружине, который ещё и тормозит о воздух.

Уравнения показаны на рисунке. Но вид системы отличается от требуемого, в том числе потому, что в нём присутствуют вторые производные. Для приведения системы в требуемый вид выполним 2 простых шага:

Первое: следует заменить переменные соответствующим образом. Теперь у нас 4 неизвестных. Далее следует преобразовать уравнение с учетом замены. Таким образом, мы имеем систему из четырёх дифференциальных уравнений первого порядка.

Настало время её записать.

Итак, мы имеем систему, параметры, интервал времени и начальные условия. Решим же эту задачу скорее.

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

Эту функцию можно располагать как в самом скрипте решения в самом его конце, так и в виде отдельного m-файла.

На выходе функция должна представлять собой вектор-столбец, который записывается перечислением компонент через точку запятой как показано на рисунке.

Теперь рассмотрим скрипт самого решения.

На этот раз запишем интервал и начальные условия в виде переменных MATLAB. Интервал, соответственно, в виде строки, а начальные условия – в виде столбца длинной 4.

Сообразно с уже разобранным ранее синтаксисом укажем функцию pendulum_np, интервал времени и начальные условия.

Перейдем теперь в окно MATLAB и посмотрим решение.

Так выглядит скрипт:

Запускаем скрипт и получаем графики:

Зачастую хочется, чтобы одну и ту же систему можно было бы решать с разными параметрами, и при этом не менять их в теле самой функции. И это можно, и даже нужно осуществлять.

На рисунке показана функция MATLAB, которая соответствует движению подвешенного на пружине шара, однако можно заметить, что эта функция теперь имеет на 5 аргументов больше.

Параметры задаются в скрипте, а при вызове функции мы обращаемся к уже известному оператору-собаке, которая превращает функцию семи переменных pendulum_n в функцию двух переменных t и X. Вот и всё.

Я вам очень рекомендую разобраться с тем, как работает оператор-собака. В хелпе он называется function-handle. Разобравшись с ним Вам будет работать в среде MATLAB ещё проще и ещё приятнее.

Вывод: не так страшно решать диффуры

Под конец стоит сказать какие вообще системы дифференциальных уравнений матлаб может решать, а может он решать системы практически любых типов.

Их можно, с одной стороны, разделить по степени жёсткости, а с другой стороны, по структуре самой системы.

Когда уравнения представляют поведение системы, которая содержит ряд быстрых и медленных реакций, то такую систему уравнения можно назвать жесткой. Для жестких задач явные численные методы работают плохо, или не работают вовсе. Примером жесткой задачи может являться протекание тока через клеточную мембрану. На самом деле, чёткого разделения между жесткими и нежёсткими системами не существует. Степень жесткости системы формально определяется через собственные значения матрицы Якоби, но давайте не будем закапываться.

Видеообзор по теме решения систем Д/У доступен по ссылке.

�� Видео

Решение систем Д/У: 1. Знакомство с функциями odeXY Скачать

Решение систем обыкновенных дифференциальных уравнений в среде MATLAB. Часть 1

В среде MATLAB можно решать системы диффуров с начальными условиями, краевые задачи, а также решать дифференциальные уравнения в частных производных с помощью инструмента PDE toolbox.

В данном обзоре речь пойдет лишь о системах дифференциальных уравнений с начальными условиями, то есть о задаче Коши. В англоязычной литературе это называется Initial Value Problem.

• каким образом записывать системы диффуров
• как задать начальные условия
• временной интервал
• какой получать результат решения для дальнего использования

Решать системы обыкновенных дифференциальных уравнений можно как в MATLAB, так и в Simulink.

В первую очередь, следует определиться, использовать для решения Matlab и его текстовый редактор, или Simulink, где те же системы дифференциальных уравнений могут быть записаны в виде функциональных блоков.

Выбор ваш должен зависеть от задачи. Если Вы, например, хотите смоделировать какой-либо объект управления, описываемый системой диффуров, то в данном случае имеет смысл использовать именно Simulink, так как Вам, впоследствии, понадобиться синтез, например, системы управления, и Simulink подойдет здесь как нельзя лучше.

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

Рассмотрим синтаксис решателей matlab.В качестве аргументов следует подать правую часть системы в виде MATLAB-функции.

На рисунке показан требуемый вид системы, когда выражены старшие производные.

Системы, чей вид отличается от требуемого, следует преобразовать к таковому.

Если функция простая, то её можно записать прямо в поле аргумента, однако, когда речь идёт о системах уравнений, имеет смысл записывать систему уравнений в виде отдельной функции, в том числе и в виде отдельного м-файла. Об этом мы поговорим чуть позже и на конкретном примере.

Также подается интервал времени, на котором будет найдено решение. Интервал задаётся строкой из двух чисел: начальной величины независимого аргумента t и его конечного значения.

Далее задаются начальные условия. Значения всех неизвестных искомых переменных в начале расчёта задаются в виде столбца соответствующей размерности.

Далее, при необходимости, задаются опции. Вот тут и раскрываются широкие возможности MATLAB по настройке решателя. Помимо управления точностью и величиной шага, имеется возможность обрабатывать данные в процессе вычисления, а также выполнять скрипты по завершению вычисления. Но ещё более полезным является опция отслеживания событий по условию, более подробно поговорим об этом дальше. Также есть другие специальные опции, которые могут быть использованы при решении определённых типов систем.

Вы могли заметить, что название функции — odeXY – это обозначение для всех решателей, которых всего 8 штук. В данном ролике мы пользоваться решателем ode45, соответствующего численному по методу Дормана-Принса 4(5). Этого решателя достаточно для подавляющего большинства задач. Остальные решатели будут подробно рассмотрены в приложении к задачам соответствующих типов позже.

Перейдем к примерам.

Рассмотрим 2 примера:

• решение дифференциального уравнения первого порядка.
• решение системы двух дифференциальных уравнений второго порядка.

В качестве уравнение первого порядка рассмотрим логистическое уравнение Ферхюльста, которое описывает динамику численности популяции. Суть уравнения такова: скорость прироста населения y пропорциональна количеству населения, однако лимитирована максимальной численностью популяции.

Забавный факт: Ферхюльст назвал это уравнение логистическим, и никто до сих пор не знает почему, ибо сам Ферхюльст об этом никому не рассказал.

Решение этого дифференциального уравнения выглядит следующим образом.

Пишем функцию в явном виде, задаём интервал расчёта и записываем начальное условие. Пару слов о записи функции подобным образом. Знак собаки в матлабе является оператором создания функции соответствующих переменных. Вы задаёте аргументы функции и саму функцию через пробел, как показано на рисунке.

Перейдем в окно MATLABа и посмотрим, как это выглядит.

Так выглядит скрипт:

Так выглядит график решения дифференциального уравнения:

В качестве примера решения системы, состоящей из двух дифференциальных уравнений второго порядка, рассмотрим шарик, подвешенный на пружине, который ещё и тормозит о воздух.

Уравнения показаны на рисунке. Но вид системы отличается от требуемого, в том числе потому, что в нём присутствуют вторые производные. Для приведения системы в требуемый вид выполним 2 простых шага:

Первое: следует заменить переменные соответствующим образом. Теперь у нас 4 неизвестных. Далее следует преобразовать уравнение с учетом замены. Таким образом, мы имеем систему из четырёх дифференциальных уравнений первого порядка.

Настало время её записать.

Итак, мы имеем систему, параметры, интервал времени и начальные условия. Решим же эту задачу скорее.

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

Эту функцию можно располагать как в самом скрипте решения в самом его конце, так и в виде отдельного m-файла.

На выходе функция должна представлять собой вектор-столбец, который записывается перечислением компонент через точку запятой как показано на рисунке.

Теперь рассмотрим скрипт самого решения.

На этот раз запишем интервал и начальные условия в виде переменных MATLAB. Интервал, соответственно, в виде строки, а начальные условия – в виде столбца длинной 4.

Сообразно с уже разобранным ранее синтаксисом укажем функцию pendulum_np, интервал времени и начальные условия.

Перейдем теперь в окно MATLAB и посмотрим решение.

Так выглядит скрипт:

Запускаем скрипт и получаем графики:

Зачастую хочется, чтобы одну и ту же систему можно было бы решать с разными параметрами, и при этом не менять их в теле самой функции. И это можно, и даже нужно осуществлять.

На рисунке показана функция MATLAB, которая соответствует движению подвешенного на пружине шара, однако можно заметить, что эта функция теперь имеет на 5 аргументов больше.

Параметры задаются в скрипте, а при вызове функции мы обращаемся к уже известному оператору-собаке, которая превращает функцию семи переменных pendulum_n в функцию двух переменных t и X. Вот и всё.

Я вам очень рекомендую разобраться с тем, как работает оператор-собака. В хелпе он называется function-handle. Разобравшись с ним Вам будет работать в среде MATLAB ещё проще и ещё приятнее.

Вывод: не так страшно решать диффуры

Под конец стоит сказать какие вообще системы дифференциальных уравнений матлаб может решать, а может он решать системы практически любых типов.

Их можно, с одной стороны, разделить по степени жёсткости, а с другой стороны, по структуре самой системы.

Когда уравнения представляют поведение системы, которая содержит ряд быстрых и медленных реакций, то такую систему уравнения можно назвать жесткой. Для жестких задач явные численные методы работают плохо, или не работают вовсе. Примером жесткой задачи может являться протекание тока через клеточную мембрану. На самом деле, чёткого разделения между жесткими и нежёсткими системами не существует. Степень жесткости системы формально определяется через собственные значения матрицы Якоби, но давайте не будем закапываться.

Видеообзор по теме решения систем Д/У доступен по ссылке.

Источники