odeint не дает мне того, что ожидал - PullRequest
Новая
       138

odeint не дает мне того, что ожидал

1 голос
/ 22 сентября 2019

Я начинаю с ODE, и у меня есть код на python, который сначала пытается решить DOE с помощью функции odeint, а затем сравнивает это решение с решением для ODE, которое вы вычислили самостоятельно (которое вы должны вставить вкод).Однако они не совсем совпадают.

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

import numpy as np
import matplotlib.pyplot as plt

# Define a function which calculates the derivative
def dy_dx(y,x):
    dydx = 2.0*x*(y-1)/(1-x**2) #
    return dydx

xs = np.linspace(-4,4,100) 

for i in range(-10,10,1):
    y0 = i
    ys = odeint(dy_dx, y0, xs) 
    ys = np.array(ys).flatten() 
    plt.plot(xs, ys);

# SECOND PART:

y_exact =  1+(y0)/(1-x**2) # manually computed solution

y_difference = ys - y_exact 
plt.subplot(2, 1, 1)
plt.plot(xs, ys, xs, y_exact, "--"); 
plt.title("Difference between exact solution and computed solution");

Итак, я добавил «часть range ()», чтобы увидеть, как она меняется в зависимости от начальных условий, но все они расположены вдоль оси x = -1

Решение, которое я нашел вручную, имеетпредел там, но это не просто строка, как вы можете видеть, запускаете ли вы вторую часть или посещаете https://www.symbolab.com/solver/ordinary-differential-equation-calculator/2x%5Cleft(y-1%5Cright)dx%2B%5Cleft(x%5E%7B2%7D-1%5Cright)dy%3D%200

Мне просто интересно, где ошибка или почему odeint дает это как результат. разницамежду odeint и моим решением

Я должен также добавить, что ODE может быть немного странным, потому что вы получаете абсолютные значения при интеграции.Это может быть связано.

1 Ответ

0 голосов
/ 22 сентября 2019

Ваше начальное условие для численного решения - y(-4)=y0, поскольку odeint принимает первую точку временного интервала в качестве начального времени.Соответственно, вам придется изменить ваше точное решение на 

y_exact =  1+((y0-1)*(1-xs[0]**2))/(1-xs**2)

, поскольку вы можете проверить с помощью инструмента по вашему выбору, например WA y'(x)=2*x*(y(x)-1)/(1-x^2), y(a)=b.

AsВы также заметили, что ваш ODE является единичным в x=-1 и x=1, так что любое решение имеет в качестве домена только один из 3 подинтервалов, созданных этими точками, а именно тот, который содержит начальное время.Кроме того, численные методы не очень хорошо работают с особенностями, поэтому вам придется ограничить область интеграции до [-4, -1-delta] для некоторых маленьких, но не слишком маленьких delta.Или, если вы действительно хотите исследовать средний интервал с начальным временем в 0, вам нужно выполнить две интеграции: одну от 0 до -1+delta и одну от 0 до 1-delta.

Если вы рассматриваете соответствующее дифференциальное уравнение без вещественных особенностей 

y'(x) = -2*x*(y-1)/(1+x^2),  y(x0) = y0

с точным решением 

y(x) = 1 + (y0-1)*(1+x0^2)/(1+x^2)

, реализованное с помощью модифицированного кода 

import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint

# Define a function which calculates the derivative
def dy_dx(y,x): return  -2.0*x*(y-1)/(1+x**2) #

# Flow function
def phi(x,x0,y0): return 1 + (y0-1)*(1+x0**2)/(1+x**2)

xs = np.linspace(-4,4,100) 

for i in range(-10,10,2):
    y0 = 1+0.1*i
    ys = odeint(dy_dx, y0, xs) 
    ys = ys.flatten() 
    plt.plot(xs, phi(xs,xs[0],y0),c='lightgray', lw=4);
    plt.plot(xs, ys, lw=1, label='$y_0=%.2f$'%y0);
plt.grid(); plt.xlim(-5.5,4.2); plt.legend(); plt.show()

, выПолучите следующий график, где вы найдете (цветное) численное решение, точно отцентрированное внутри (подкладка серого) точного решения.enter image description here

...