Блог

Итак, мне было интересно: не должна ли сумма интегралов некоторых функций быть эквивалентна интегралу от суммы функций?

Здесь я интегрирую три произвольные функции с помощью scipy ’s solve_ivp :

 import numpy as np
from scipy.integrate import solve_ivp

def fun1(t, y): return 0.5 * y
def fun2(t, y): return 0.05 * y**2
def fun3(t, y): return y   5

sol1 = solve_ivp(fun1, [0, 5], [2], rtol=1e-10, atol=1e-10)
sol2 = solve_ivp(fun2, [0, 5], [2], rtol=1e-10, atol=1e-10)
sol3 = solve_ivp(fun3, [0, 5], [2], rtol=1e-10, atol=1e-10)

print(f'{sol1.y[0][-1]}   {sol2.y[0][-1]}   {sol3.y[0][-1]} = {sol1.y[0][-1] sol2.y[0][-1] sol3.y[0][-1]}')
  

это дает:

 24.36498792283038   4.0000000000803775   1033.8921138100272 = 1062.257101732938
  

Но, с другой стороны:

 import numpy as np
from scipy.integrate import solve_ivp

def fun(t, y): return (0.5 * y)   (0.05 * y**2)   (y   5)

sol = solve_ivp(fun, [0, 5], [2], rtol=1e-10, atol=1e-10)

print(f'{sol.y[0][-1]}')
  

дает:

 233176411972824.97
  

Может ли кто-нибудь из вас, умных людей, показать мне мою ошибку в мышлении? Моя математика неверна или я должен реализовать эту задачу по-другому в Python? Большое вам спасибо за вашу помощь!

1. Начальные значения несовместимы. Вы также должны суммировать начальные значения.
2. Если вы суммируете решения

 dy/dt = -y 
dy/dt = y 
y(0)=1
  

вы не получаете решение

 dy/dt = 0 
  

вы получаете

 cos(t)   y^t
  

Но

 dy/dt = t
dy/dt = t 
  

вы действительно получаете решение

 dy/dt = 2t 
  

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

 import numpy as np
from scipy.integrate import solve_ivp

def fun1(t, y): return 0.5 * t
def fun2(t, y): return 0.05 * t**2
def fun3(t, y): return t   5

sol1 = solve_ivp(fun1, [0, 5], [2], rtol=1e-10, atol=1e-10)
sol2 = solve_ivp(fun2, [0, 5], [2], rtol=1e-10, atol=1e-10)
sol3 = solve_ivp(fun3, [0, 5], [2], rtol=1e-10, atol=1e-10)

print(f'{sol1.y[0][-1]}   {sol2.y[0][-1]}   {sol3.y[0][-1]} = {sol1.y[0][-1] sol2.y[0][-1] sol3.y[0][-1]}')


def fun(t, y): return (0.5 * t)   (0.05 * t**2)   (t   5)

sol = solve_ivp(fun, [0, 5], [6], rtol=1e-10, atol=1e-10)

print(f'{sol.y[0][-1]}')
  

Обновить

Как правило, вы решаете dy/dt=f(y,t) . Случай dy/dt=g(t) отличается от dy/dt=h(y) .

Если у вас нет y в правой части, вы можете просто суммировать уравнения

 dy1/dt=g1(t) 
dy2/dt=g2(t)
=>
d(y1   y2)/dt=g1(t)   g2(t)
  

сделайте замену y = dy1 dy2 и получите

 y/dt=g1(t)   g2(t)
  

И все идет так, как вы ожидаете. Вы можете применить solve_ivp ко всей функции или к ее частям.

Если у вас есть y в правой части

 dy1/dt=h1(y1)
dy2/dt=h2(y2)
  

и сумма, которую вы получаете

 d(y1   y2)/dt=h1(y1)   h2(y2)
  

и подстановка y = y1 y2 никуда вас не приведет. Таким образом, вы не можете разделить правую часть для такого случая на части, чтобы получить весь результат путем отдельной интеграции.