python 经典R-K(龙格库塔)

经典R-K(龙格库塔)

import sympy
x=sympy.symbols("x")
y=sympy.symbols("y")
fy=input("f(x,y)=")
a=float(input("a="))
b=float(input("b="))
n=int(input("n="))
α=float(input("α="))
fy=sympy.simplify(fy)
X=[]
Y=[]
h=(b-a)/n
X.append(a)
Y.append(α)
for i in range(n):
    x1=X[i]+h
    X.append(x1)
    K1=fy.evalf(subs={x:X[i],y:Y[i]})       #对x,y赋值
    K2=fy.evalf(subs={x:(X[i]+h/2),y:(Y[i]+(h/2)*K1)})
    K3=fy.evalf(subs={x:(X[i]+h/2),y:(Y[i]+(h/2)*K2)})
    K4=fy.evalf(subs={x:X[i+1],y:(Y[i]+h*K3)})
    y1=Y[i]+h*(1/6)*(K1+2*K2+2*K3+K4)
    Y.append(y1)
print('x:',X)
print('y:',Y)

输入案例：