как решить дифференциальное уравнение в Python

Question

y''' + yy" + (1 - y'^2)= 0 y(0)=0, y'(0)=0, y'(+∞)=1 

(The +∞ может быть заменен на 4) это falkenr-скан equaiton. Я хочу получить решение numreiacl от 0 до 4.

На самом деле, я прочитал книгу «Численные методы в технике с Python 3» и использовал методы в книге, но я не могу получить anwser, может быть, там некоторые ошибки в книге. здесь код

""" 
solve the differential equation 
y''' + yy" + (1 - y'^2)= 0 y(0) = 0 y'(0) = 0,y'(+inf) = 0 
""" 
import numpy as np 
from scipy.integrate import odeint 
from printSoln import * 
start = time.clock() 
def F(x,y): # First-order differential equations 
    y = np.zeros(3) 
    F0 = y[0] 
    F1 = y[2] 
    F2 = y[2] 
    F3 = -y[0]*y[2] - (1 - y[1]*y[1]) 
    return [F1,F2,F3] 

def init(u): 
    return np.array([0,0,u]) 
X = np.linspace(0, 5, 50) 
u = 0 
for i in range(1): 
    u += 1 
    Y = odeint(F, init(u), X) 
    print(Y) 
    if abs(Y[len(Y)-1,1] - 1) < 0.001: 
     print(i,'times iterations') 
     printSoln(X,Y,freq=2) 
     print("y'(inf)) = ", Y[len(Y)-1,1]) 
     print('y"(0) =',u) 
     break 
end = time.clock() 

print('duration =',end - start) 

Answer 1

код вы показываете, как предполагается реализовать метод съемки, чтобы решить краевую задачу, сведя ее к начальной задачи (https://en.wikipedia.org/wiki/Shooting_method). Однако в коде есть много ошибок.

1. Читая документацию (https://docs.scipy.org/doc/scipy-0.18.1/reference/generated/scipy.integrate.odeint.html) для odeint функции, которая является ядром решения, можно обнаружить, что первый аргумент odeint является func(y, t0, ...). Поэтому мы должны заменить линию def F(x,y): на def F(y,x):

2. Функция F(y, x) должна вычислять производную от y по x. Переписывая исходный 3-й порядок нелинейного ОДУ как систему из трех ОДЕ 1-го порядка, можно найти, что: (а) строка y = np.zeros(3) не имеет смысла, поскольку она заставляет все производные быть нулевыми; (b) оно должно быть F1 = y[1] вместо F1 = y[2] и (c) строка F0 = y[0] не является обязательной и может быть удалена.

3. Мы уже выяснили, что код должен реализовывать метод съемки. Это означает, что на левой границе мы должны установить условия y(0) = 0, y'(0) = 0, которые мы имеем из утверждения проблемы. Однако для решения проблемы начального значения нам нужно еще одно условие для y''(0). Идея состоит в том, что мы решим нашу систему ОДУ с различными условиями вида y''(0) = u, пока для некоторого значения u решение удовлетворяет граничному условию y'(4) = 1 на правой границе с заданным допуском. Следовательно, для целей экспериментов заменим строку for i in range(1): на for i in range(5): и напечатаем значение y'(4) на print Y[-1,1]. Выход будет:

   -1.26999326625 
   19.263464565 
   73.5661968483 
   175.047093183 
   340.424666137 
   

Можно видеть, что если приращение u = 1: (а) y'(4) монотонно возрастать и (б) при и = 1 y'(4)=-1.26999326625, при и = 2 y'(4)=19.263464565. Поэтому решение, которое удовлетворяет граничному условию y'(4)=1, может быть найдено с 1<u<2. Таким образом, мы должны уменьшить приращение u, чтобы найти решение с более высоким допуском.

Финлей, мы можем переписать код следующим образом:

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

start = time.clock() 
def F(y, x): # First-order differential equations 
    F1 = y[1] 
    F2 = y[2] 
    F3 = -y[0]*y[2] - (1 - y[1]*y[1]) 
    return [F1,F2,F3] 

def init(u): 
    return np.array([0,0,u]) 

X = np.linspace(0, 4, 50) 
niter = 100 
u = 0 
tol = 0.1 
ustep = 0.01 

for i in range(niter): 
    u += ustep 
    Y = odeint(F, init(u), X) 
    print Y[-1,1] 
    if abs(Y[len(Y)-1,1] - 1) < tol: 
     print(i,'times iterations') 
     print("y'(inf)) = ", Y[len(Y)-1,1]) 
     print('y"(0) =',u) 
     break 

end = time.clock() 
print('duration =',end - start) 

plt.plot(X, Y[:,0], label='y') 
plt.plot(X, Y[:,1], label='y\'') 
plt.legend() 
plt.show() 

Численные исследования с ustep=0.01 показывает, что два Diferent решения возможны. Один на 0<u<1 и другой на 1<u<2. Эти решения показаны ниже.

Answer 2

# biscetion 
""" root = bisection(f,x1,x2,switch=0,tol=1.0e-9). 
Finds a root of f(x) = 0 by bisection. 
The root must be bracketed in (x1,x2). 
Setting switch = 1 returns root = None if 
f(x) increases upon bisection. 
""" 
import math 
from numpy import sign 
def bisection(f,x1,x2,switch=1,tol=1.0e-9): 
    f1 = f(x1) 
    if f1 == 0.0: return x1 
    f2 = f(x2) 
    if f2 == 0.0: return x2 
    if sign(f1) == sign(f2): 
     print('Root is not bracketed') 
    n = int(math.ceil(math.log(abs(x2 - x1)/tol)/math.log(2.0))) 
    for i in range(n): 
     x3 = 0.5*(x1 + x2); f3 = f(x3) 
     if (switch == 1) and (abs(f3) > abs(f1)) \ 
      and (abs(f3) > abs(f2)): 
      return None 

     if f3 == 0.0: return x3 
     if sign(f2)!= sign(f3): x1 = x3; f1 = f3 
     else: x2 = x3; f2 = f3 
    return (x1 + x2)/2.0 
""""//////////////////////////////////////////////////////////////////""" 
# run_kut4 
""" X,Y = integrate(F,x,y,xStop,h). 
    4th-order Runge-Kutta method for solving the 
    initail value problem {y}' = {F(x,{y})},where 
    {y} = {y[0],y[1],...,y[n-1]}. 
    x,y = initial conditions 
    xStop = terminal value of x 
    h  =increment of x used in integration 
    F  =user - supplied function that returns the 
      array F(x,y) = {y'[0],y'[1],...,y'[n-1]}. 
""" 
import numpy as np 
def integrate(F,x,y,xStop,h): 

    def run_kut4(F,x,y,h): 
     K0 = h*F(x,y) 
     K1 = h*F(x + h/2.0, y + K0/2.0) 
     K2 = h*F(x + h/2.0, y + K1/2.0) 
     K3 = h*F(x + h, y + K2) 
     return (K0 + 2.0*K1 + 2.0*K2 + K3)/6.0 
    X = [] 
    Y = [] 
    X.append(x) 
    Y.append(y) 
    while x < xStop: 
     h = min(h,xStop - x) 
     y = y + run_kut4(F,x,y,h) 
     x = x + h 
     X.append(x) 
     Y.append(y) 
    return np.array(X), np.array(Y) 
""""//////////////////////////////////////////////////////////////////"""  
## printsSoln 
""" printSoln(X,Y,freq). 
    Prints X and Y returner from the differential 
    equation solves using printout frequency 'freq'. 
     freq = n prints every nth step. 
     freq = 0 prints inital and final values only. 
""" 
def printSoln(X,Y,freq): 

    def printHead(n): 
     print('\n{:>13}'.format('x'),end=" ") 
     for i in range(n): 
      print('{}{}{}{}'.format(' '*(9),'y[',i,']'),end=" ") 
     print() 
    def printLine(x,y,n): 
     print("{:13.4}".format(x),end=" ") 
     for i in range(n): 
      print("{:13.4f}".format(y[i]),end=" ") 
     print() 
    m = len(Y) 
    try: n = len(Y[0]) 
    except TypeError: n = 1 
    if freq == 0: freq = m 
    printHead(n) 
    for i in range(0,m,freq): 
     printLine(X[i],Y[i],n) 
    if i != m - 1: printLine(X[m - 1],Y[m - 1],n) 
""""//////////////////////////////////////////////////////////////////""" 

""" 
solve the differential equation 
y''' + yy" + B(1 - y'^2)= 0 y(0) = 0 y'(0) = 0,y'(+inf) = 0 
B = [0,0.2,0.4,0.6,0.8,1.0,-0.1,-0.15] 
""" 
import matplotlib.pyplot as plt 
import time 
start = time.clock() 
def initCond(u): #Init. values of [y.y']; use 'u' if unkown 
    return np.array([0.0,0.0,u]) 

def r(u): # Boundart condition resdiual 
    X,Y = integrate(F,xStart,initCond(u),xStop,h) 
    y = Y[len(Y) - 1] 
    r = y[1] - 1.0 
    return r 

def F(x,y): # Third-order differential equations 
    F = np.zeros(3) 
    F[0] = y[1] 
    F[1] = y[2] 
    F[2] = -y[0]*y[2] - (1 - y[1]*y[1]) 
    return F 

xStart = 0.0 # Start of integration 
xStop = 5 # End of integraiton 
u1 = 1.0  # 1st trial value of unknown init. cond. 
u2 = 2.0  # 2nd trial value of unknown init. cond. 
h = 0.1  # Step size 
freq = 2  # Printout frequency 


u = bisection(r,u1,u2,tol = 1.0e-9) #Comput the correct initial condition 
X,Y = integrate(F,xStart,initCond(u),xStop,h) 
print(initCond(u)) 
printSoln(X,Y,freq) 
end = time.clock() 
print('duration =',end - start) 
plt.plot(X, Y[:,0], label='y') 
plt.plot(X, Y[:,1], label='y\'') 
plt.legend() 
plt.show()