Хабрахабр

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

Введение:

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

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

Что касается Python, то в публикациях по численным методам, например [2,3], данных по применение Рунге — Кутты крайне мало, а по его модификации — методу Рунге-Кутта-Фельберга вообще нет. Возникает необходимость использовать численные методы, наиболее известным из которых является метод Рунге — Кутты [1].

Вторая функция ode из этого модуля реализует несколько методов, в том числе и упомянутый пятиранговый метод Рунге-Кутта-Фельберга, но, вследствие универсальности, имеет ограниченное быстродействие. В настоящее время, благодаря простому интерфейсу, наибольшее распространение в Python имеет функцию odeint из модуля scipy.integrate.

В публикации так же приведены решения по краевым задачам для систем дифференциальных уравнений (СДУ).
Целью настоящей публикации является сравнительный анализ перечисленных средств численного решения систем дифференциальных уравнений с модифицированным автором под Python методом Рунге-Кутта-Фельберга.

Краткие теоретические и фактические данные по рассматриваемым методам и программным средствам для численного решения СДУ

Задача Коши

Для одного дифференциального уравнения n – го порядка, задача Коши состоит в нахождении функции, удовлетворяющей равенству:

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

Перед решением эта задача должна быть переписана в виде следующей СДУ

(1)

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

Модуль scipy.integrate

Функция ode() более универсальная, а функция odeint() (ODE integrator) имеет более простой интерфейс и хорошо решает большинство задач. Модуль имеет две функции ode() и odeint(), предназначенные для решения систем обыкновенных дифференциальных уравнений (ОДУ) первого порядка с начальными условиями в одной точке (задача Коши).

Функция odeint()

Она имеет следующий формат odeint(func, y0, t[,args=(), ...]) Аргумент func – это имя Python функции двух переменных, первой из которых является список y=[y1,y2,...,yn], а второй – имя независимой переменной. Функция odeint() имеет три обязательных аргумента и много опций.

Фактически функция func(y,t) реализует вычисление правых частей системы (1). Функция func должна возвращать список из n значений функций при заданном значении независимого аргумента t.

Второй аргумент y0 функции odeint() является массивом (или списком) начальных значений при t=t0.

При этом первый элемент этого массива рассматривается как t0. Третий аргумент является массивом моментов времени, в которые вы хотите получить решение задачи.

Функция odeint() имеет много опций, управляющих ее работой. Функция odeint() возвращает массив размера len(t) x len(y0). Опции rtol (относительная погрешность) и atol (абсолютная погрешность) определяют погрешность вычислений ei для каждого значения yi по формуле

По умолчанию Они могут быть векторами или скалярами.

Функция ode()

Она создает объект ОДУ (тип scipy.integrate._ode.ode). Вторая функция модуля scipy.integrate, которая предназначена для решения дифференциальных уравнений и систем, называется ode(). Аналогично функции odeint(), функция ode(func) предполагает приведение задачи к системе дифференциальных уравнений вида (1) и использовании ее функции правых частей. Имея ссылку на такой объект, для решения дифференциальных уравнений следует использовать его методы.

Например, следующая последовательность инструкций создает объект ODE, представляющий задачу Коши. Отличие только в том, что функция правых частей func(t,y) первым аргументом принимает независимую переменную, а вторым – список значений искомых функций.

Метод Рунге—Кутта

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

При численном решении задачи Коши

(2)

(3)

При численном решении задачи (2),(3) будем использовать равномерную, для простоты, сетку по переменной t с шагом т > 0. по известному решению в точке t =0 необходимо найти из уравнения (3) решение при других t.

Метод сходится в точке если при . Приближенное решение задачи (2), (3) в точке обозначим . Простейшая разностная схема для приближенного решения задачи (2),(3) есть Метод имеет р-й порядок точности, если , р > 0 при .

(4)

Симметричная схема в (4) имеет второй порядок аппроксимации. При имеем явный метод и в этом случае разностная схема аппроксимирует уравнение (2) с первым порядком. Эта схема относится к классу неявных — для определения приближенного решения на новом слое необходимо решать нелинейную задачу.

На этапе предиктора (предсказания) используется явная схема Явные схемы второго и более высокого порядка аппроксимации удобно строить, ориентируясь на метод предиктор-корректор.

(5)

а на этапе корректора (уточнения) — схема

Этот метод записывается в общем виде: В одношаговых методах Рунге—Кутта идеи предиктора-корректора реализуются наиболее полно.

(6),

где

Если при имеем явный метод Рунге—Кутта. Формула (6) основана на s вычислениях функции f и называется s-стадийной. Если при j>1 и то определяется неявно из уравнения:

(7)

Параметры определяют вариант метода Рунге—Кутта. О таком методе Рунге—Кутта говорят как о диагонально-неявном. Используется следующее представление метода (таблица Бутчера)

Одним из наиболее распространенных является явный метод Рунге—Кутта четвертого порядка

(8)

Метод Рунге—Кутта— Фельберга

Привожу значение расчётных коэффициентов метода

(9)

С учётом(9) общее решение имеет вид:

(10)

Это решение обеспечивает пятый порядок точности, остаётся его адаптировать к Python.

Вычислительный эксперимент по определению абсолютной погрешности численного решения нелинейного дифференциального уравнения с использованием обеих функций def odein(),def oden() модуля scipy.integrate и адаптированного к Python методов Рунге—Кутта и Рунге—Кутта— Фельберга

Листинг программы

from numpy import*
import matplotlib.pyplot as plt
from scipy.integrate import *
def odein(): #dy1/dt=y2 #dy2/dt=y1**2+1: def f(y,t): return y**2+1 t =arange(0,1,0.01) y0 =0.0 y=odeint(f, y0,t) y = array(y).flatten() return y,t
def oden(): f = lambda t, y: y**2+1 ODE=ode(f) ODE.set_integrator('dopri5') ODE.set_initial_value(0, 0) t=arange(0,1,0.01) z=[] t=arange(0,1,0.01) for i in arange(0,1,0.01): ODE.integrate(i) q=ODE.y z.append(q[0]) return z,t def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau): if z==1: k0 =tau* f(t,y) k1 =tau* f(t+tau/2.,y+k0/2.) k2 =tau* f(t+tau/2.,y+k1/2.) k3 =tau* f(t+tau, y + k2) return (k0 + 2.*k1 + 2.*k2 + k3) / 6. elif z==0: k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = [] y= [] t.append(to) y.append(yo) while to < tEnd: tau = min(tau, tEnd - to) yo = yo + increment(f, to, yo, tau) to = to + tau t.append(to) y.append(yo) return array(t), array(y)
def f(t, y): f = zeros([1]) f[0] = y[0]**2+1 return f
to = 0.
tEnd = 1
yo = array([0.])
tau = 0.01
z=1
t, yn = rungeKutta(f, to, yo, tEnd, tau)
y1n=[i[0] for i in yn]
plt.figure()
plt.title("Абсолютная погрешность численного решения(т.р.- u(t)=tan(t)) ДУ\n\
du/dt=u**2+1 c u(0)=0 при t>0")
plt.plot(t,abs(array(y1n)-array(tan(t))),label='Метод Рунге—Кутта \n\
четвертого порядка - расчёт по алгоритму')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
z=0
t, ym = rungeKutta(f, to, yo, tEnd, tau)
y1m=[i[0] for i in ym]
plt.figure()
plt.title("Абсолютная погрешность численного решения(т.р.- u(t)=tan(t)) ДУ\n\
du/dt=u**2+1 c u(0)=0 при t>0")
plt.plot(t,abs(array(y1m)-array(tan(t))),label='Метод Рунге—Кутта— Фельберга \n\
пятого порядка - расчёт по алгоритму')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
plt.figure()
plt.title("Абсолютная погрешность численного решения (т.р.- u(t)=tan(t)) ДУ\n\
du/dt=u**2+1 c u(0)=0 при t>0")
y,t=odein()
plt.plot(t,abs(array(tan(t))-array(y)),label='Функция odein')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
plt.figure()
plt.title("Абсолютная погрешность численного решения (т.р.- u(t)=tan(t)) ДУ\n\
du/dt=u**2+1 c u(0)=0 при t>0")
z,t=oden()
plt.plot(t,abs(tan(t)-z),label='Функция ode метод Рунге—Кутта— Фельберга \n\
пятого порядка')
plt.xlabel('Время')
plt.ylabel('Абсолютная погрешность.')
plt.legend(loc='best')
plt.grid(True)
plt.show()

Получим:

Вывод:

Необходимо провести исследование быстродействия. Адаптированные к Python методы Рунге—Кутта и Рунге—Кутта— Фельберга имеют меньшую абсолютную, чем решение с применением функции odeint, но большую, чем с использованием функции edu.

Численный эксперимент по сравнению быстродействия численного решения СДУ при использовании функции ode с атрибутом dopri5 (метод Рунге – Кутты 5 порядка) и с использованием адаптированного к Python метода Рунге—Кутта— Фельберга

Чтобы не повторять источник, приведу постановку и решение модельной задачи из [2]. Сравнительный анализ проведём на примере модельной задачи, приведенной в [2].

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

где – радиус вектор движущегося тела, – вектор скорости тела, – коэффициент сопротивления, вектор силы веса тела массы m, g – ускорение свободного падения.

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

Положим . К системе следует добавить начальные условия: (h начальная высота), . Тогда соответствующая система ОДУ 1 – го порядка примет вид:

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

Листинг программы

import numpy as np
import matplotlib.pyplot as plt
import time
start = time.time()
from scipy.integrate import ode
ts = [ ]
ys = [ ]
FlightTime, Distance, Height =0,0,0
y4old=0
def fout(t, y):# обработчик шага global FlightTime, Distance, Height,y4old ts.append(t) ys.append(list(y.copy())) y1, y2, y3, y4 = y if y4*y4old<=0: # достигнута точка максимума Height=y3 if y4<0 and y3<=0.0: # тело достигло поверхности FlightTime=t Distance=y1 return -1 y4old=y4
# функция правых частей системы ОДУ
def f(t, y, k): # имеется дополнительный аргумент k g=9.81 y1, y2, y3, y4 = y return [y2,-k*y2*np.sqrt(y2**2+y4**2), y4,-k*y4*np.sqrt(y2**2+y4**2)-g] tmax=1.41 # максимально допустимый момент времени
alph=np.pi/4 # угол бросания тела
v0=10.0 # начальная скорость
K=[0.1,0.2,0.3,0.5] # анализируемые коэффициенты сопротивления
y0,t0=[0, v0*np.cos(alph), 0, v0*np.sin(alph)], 0 # начальные условия ODE=ode(f)
ODE.set_integrator('dopri5', max_step=0.01)
ODE.set_solout(fout)
fig, ax = plt.subplots()
fig.set_facecolor('white')
for k in K: # перебор значений коэффициента сопротивления ts, ys = [ ],[ ] ODE.set_initial_value(y0, t0) # задание начальных значений ODE.set_f_params(k) # передача дополнительного аргумента k # в функцию f(t,y,k) правых частей системы ОДУ ODE.integrate(tmax) # решение ОДУ print('Flight time = %.4f Distance = %.4f Height =%.4f '% (FlightTime,Distance,Height)) Y=np.array(ys) plt.plot(Y[:,0],Y[:,2],linewidth=3,label='k=%.1f'% k)
stop = time.time()
plt.title("Результаты численного решения системы четырёх ОДУ \n с использованием функци ode с атрибутом dopri5 ")
print ("Время на модельную задачу: %f"%(stop-start))
plt.grid(True)
plt.xlim(0,8)
plt.ylim(-0.1,2)
plt.legend(loc='best') plt.show()

Получим:

2316 Distance = 5. Flight time = 1. 8542
Flight time = 1. 9829 Height =1. 3830 Height =1. 1016 Distance = 4. 0197 Distance = 3. 5088
Flight time = 1. 2912
Flight time = 0. 5265 Height =1. 5842 Height =1. 9068 Distance = 2. 454787 0240
Время на модельную задачу: 0.

Для реализации средствами Python численного решения СДУ без использования специальных модулей, мною была предложена и исследована следующая функция:

def increment(f, t, y, tau
k1=tau*f(t,y)
k2=tau*f(t+(1/4)*tau,y+(1/4)*k1)
k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2)
k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3)
k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4)
k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5)
return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6

Остальные особенности программы можно посмотреть в следующем листинге: Функция increment(f, t, y, tau) обеспечивает пятый порядок численного метода решения.

Листинг программы

from numpy import*
import matplotlib.pyplot as plt
import time
start = time.time()
def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau):# поиск приближённого решения методом Рунге—Кутта—Фельберга. k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = []#подготовка пустого списка t y= []#подготовка пустого списка y t.append(to)#внесение в список t начального значения to y.append(yo)#внесение в список y начального значения yo while to < tEnd:#внесение результатов расчёта в массивы t,y tau = min(tau, tEnd - to)#определение минимального шага tau yo = yo + increment(f, to, yo, tau) # расчёт значения в точке t0,y0 для задачи Коши to = to + tau # приращение времени t.append(to) # заполнение массива t y.append(yo) # заполнение массива y return array(t), array(y)
def f(t, y): # функция правых частей системы ОДУ f = zeros([4]) f[0]=y[1] f[1]=-k*y[1]*sqrt(y[1]**2+y[3]**2) f[2]=y[3] f[3]=-k*y[3]*sqrt(y[1]**2+y[3]**2) -g if y[3]<0 and y[2]<=0.0: # тело достигло поверхности return -1 return f
to = 0# начальный момент отсчёта времени
tEnd = 1.41# конечный момент отсчёта времени
alph=pi/4# угол бросания тела
v0=10.0 # начальная скорость
K=[0.1,0.2,0.3,0.5]# анализируемые коэффициенты сопротивления среды
g=9.81 yo = array([0.,v0*cos(alph),0.,v0*sin(alph)]) # начальные условия tau =0.01# шаг
for i in K: # перебор значений коэффициента сопротивления среды k=i t, y = rungeKutta(f, to, yo, tEnd, tau) y1=array([i[0] for i in y]) # извлечение переменных из массива y y3=array([i[2] for i in y]) # визуализация сопротивление среды "к" и расчёт и визуализация s,h,t plt.plot(y1,y3,linewidth=2,label='k=%.1f h=%.3f s=%.2f t=%s' % (k,max(y3),max(y1),round(t[list(y1).index(max(y1))],3)))
stop = time.time()
plt.title("Результаты численного решения системы четырёх ОДУ \n с использованием адаптированного для Python\n метода Рунге—Кутта—Фельберга ")
print ("Время на модельную задачу: %f"%(stop-start))
plt.xlabel('Высота h')
plt.ylabel('Расстояние s')
plt.legend(loc='best') plt.xlim(0,8)
plt.ylim(-0.1,2) plt.grid(True) plt.show()

Получим:

259927 Время на модельную задачу: 0.

Вывод

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

Решение краевой задачи с поточно разделёнными краевыми условиями

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

(11)

Для решения задачи (11) используем следующий алгоритм:

Решаем первые три неоднородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:
1.

Решаем первые три однородные уравнения системы (11) с начальными условиями

Введем обозначение для решения задачи Коши:
2.

Решаем первые три однородные уравнения системы (11) с начальными условиями 3.

Введем обозначение для решения задачи Коши:

Общее решение краевой задачи (11) при помощи решений задач Коши записывается в виде линейной комбинации решений:

где p2, p3 — некоторые неизвестные параметры. 4.

Для определения параметров p2, p3, используем краевые условия последних двух уравнений (11), то есть условия при x = b. 5. Подставляя, получим систему линейных уравнений относительно неизвестных p2, p3:
(12)
Решая (12), получим соотношения для p2, p3.

По приведенному алгоритму с применением метода Рунге—Кутта—Фельберга получим следующую программу:

Листинг программы

# Метод пристрелки
from numpy import*
import matplotlib.pyplot as plt
import matplotlib.font_manager as fm,os
import matplotlib.patches as mpatches
import matplotlib.lines as mlines
from scipy.integrate import odeint
from scipy import linalg
import time
start = time.time()
c1 = 1.0
c2 = 0.8
c3 = 0.5
a =0.0
b = 1.0
nn =100
initial_state_0 =array( [a, c1, 0.0, 0.0])
initial_state_I =array( [a, 0.0, 1.0, 0.0])
initial_state_II =array( [a, 0.0, 0.0, 1.0])
to = a
tEnd =b
N = int(nn)
tau=(b-a)/N
def rungeKutta(f, to, yo, tEnd, tau): def increment(f, t, y, tau): k1=tau*f(t,y) k2=tau*f(t+(1/4)*tau,y+(1/4)*k1) k3 =tau *f(t+(3/8)*tau,y+(3/32)*k1+(9/32)*k2) k4=tau*f(t+(12/13)*tau,y+(1932/2197)*k1-(7200/2197)*k2+(7296/2197)*k3) k5=tau*f(t+tau,y+(439/216)*k1-8*k2+(3680/513)*k3 -(845/4104)*k4) k6=tau*f(t+(1/2)*tau,y-(8/27)*k1+2*k2-(3544/2565)*k3 +(1859/4104)*k4-(11/40)*k5) return (16/135)*k1+(6656/12825)*k3+(28561/56430)*k4-(9/50)*k5+(2/55)*k6 t = [] y= [] t.append(to) y.append(yo) while to < tEnd: tau = min(tau, tEnd - to) yo = yo + increment(f, to, yo, tau) to = to + tau t.append(to) y.append(yo) return array(t), array(y)
def f(t, y): global theta f = zeros([4]) f[0] = 1 f[1] = -y [1]-y[2] +theta* sin(y[0]) f[2] = -y[2]+y[3] f[3] = -y[2] return f
# Решение НЕОДНОРОДНОЙ системы -- theta = 1
theta = 1.0
yo =initial_state_0
t, y = rungeKutta(f, to, yo, tEnd, tau)
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Извлечение требуемых для решения задачи значений # Y20 = Y2(b), Y30 = Y3(b)
Y20 = y2[N-1]
Y30 = y3[N-1]
# Решение ОДНОРОДНОЙ системы -- theta = 0, задача I
theta = 0.0
yo= initial_state_I
t, y = rungeKutta(f, to, yo, tEnd, tau)
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Извлечение требуемых для решения задачи значений # Y21= Y2(b), Y31 = Y3(b)
Y21= y2[N-1]
Y31 = y3[N-1]
# Решение ОДНОРОДНОЙ системы -- theta = 0, задача II
theta = 0.0
yo =initial_state_II
t, y = rungeKutta(f, to, yo, tEnd, tau)
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Извлечение требуемых для решения задачи значений # Y211= Y2(b), Y311 = Y3(b)
Y211= y2[N-1]
Y311 = y3[N-1]
# Формирование системы линейных
# АЛГЕБРАИЧЕСКИХ уравнений для определния p2, p3
b1 = c2 - Y20
b2 = c3 - Y30
A = array([[Y21, Y211], [Y31, Y311]])
bb = array([[b1], [b2]]) # Решение системы
p2, p3 = linalg.solve(A, bb)
# Окончательное решение краевой
# НЕОДНОРОДНОЙ задачи, theta = 1
theta = 1.0
yo = array([a, c1, p2, p3])
t, y = rungeKutta(f, to, yo, tEnd, tau)
y0=[i[0] for i in y]
y1=[i[1] for i in y]
y2=[i[2] for i in y]
y3=[i[3] for i in y]
# Проверка
print('y0[0]=', y0[0])
print('y1[0]=', y1[0])
print('y2[0]=', y2[0])
print('y3[0]=', y3[0])
print('y0[N-1]=', y0[N-1])
print('y1[N-1]=', y1[N-1])
print('y2[N-1]=', y2[N-1])
print('y3[N-1]=', y3[N-1])
j = N
xx = y0[:j]
yy1 = y1[:j]
yy2 = y2[:j]
yy3 = y3[:j]
stop = time.time()
print ("Время на модельную задачу: %f"%(stop-start))
plt.subplot(2, 1, 1)
plt.plot([a], [c1], 'ro')
plt.plot([b], [c2], 'go')
plt.plot([b], [c3], 'bo')
plt.plot(xx, yy1, color='r') #Построение графика
plt.plot(xx, yy2, color='g') #Построение графика
plt.plot(xx, yy3, color='b') #Построение графика
plt.xlabel(r'$x$') #Метка по оси x в формате TeX
plt.ylabel(r'$y_k(x)$') #Метка по оси y в формате TeX
plt.title(r'Метод пристрелки ', color='blue')
plt.grid(True) #Сетка
patch_y1 = mpatches.Patch(color='red', label='$y_1$')
patch_y2 = mpatches.Patch(color='green', label='$y_2$')
patch_y3 = mpatches.Patch(color='blue', label='$y_3$')
plt.legend(handles=[patch_y1, patch_y2, patch_y3])
ymin, ymax = plt.ylim()
xmin, xmax = plt.xlim()
plt.subplot(2, 1, 2)
font =
plt.text(0.2, 0.8, r'$\frac{dy_1}{dx}= - y_1 - y_2 + \sin(x),$', fontdict=font)
plt.text(0.2, 0.6,r'$\frac{dy_2}{dx}= - y_1 + y_3,$', fontdict=font)
plt.text(0.2, 0.4, r'$\frac{dy_3}{dx}= - y_2 - y_2,$', fontdict=font)
plt.text(0.2, 0.2, r'$y_1(a)=c_1, ' r'\quad y_2(b)=c_2, \quad y_3(b)=c_3.$', fontdict=font)
plt.subplots_adjust(left=0.15)
plt.show()

Получим:

0
y1[0]= 1. y0[0]= 0. 7156448588231397
y3[0]= 1. 0
y2[0]= 0. 9900000000000007
y1[N-1]= 0. 324566562303714
y0[N-1]= 0. 8
y3[N-1]= 0. 1747719838716767
y2[N-1]= 0. 070878 5000000000000001
Время на модельную задачу: 0.

Вывод

Разработанная мною программа отличается от приведенной в [3] меньшей погрешностью, что подтверждает приведенный в начале статьи сравнительный анализ функции odeint с реализованным на Python метода Рунге—Кутта—Фельберга.

Ссылки:

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

Введение в научный Python. 2.

Н.М. 3. Ширяева Python 3. Полякова, Е.В. Ростов-на-Дону 2017. Создание графического интерфейса пользователя (на примере решения методом пристрелки краевой задачи для линейных обыкновенных дифференциальных уравнений).

Теги
Показать больше

Похожие статьи

Добавить комментарий

Ваш e-mail не будет опубликован. Обязательные поля помечены *

Кнопка «Наверх»
Закрыть