python——求解微分方程

1.用预报-校正法解下面常微分方程的初值问题。

 2.用r-k方法求解初值问题：
 

---

import math as mh
import matplotlib.pyplot as plt
def g(x):
    ni=1/((x**2)+1)
    return ni
def no(jk1,h2):#question one with prediction correction method
    global y
    def f(x,y):
        xv=-2*x*(y**2)
        return xv
    x1=jk1;x2=h2;h=0.01;y0=1
    t=(x2-x1)/h;hi=int(t)
    for i in range(hi):
        n1=f(x1+i*h,y0);
        y1=y0+h*n1;
        n2=f(x1+(i+1)*h,y1);
        y=y0+h/2*(n1+n2)
        y0=y
    return y
x0=0#x的初值点
x1=2#想求的X值
r=1/5;y=0;no(x0,x1)
print('问题1:(保留6位小数)\n求解X={}时,y值\n预报校正法求解：{:.6};步长为0.01'.format(x1,y))
print('准确求解微分方程：{:.6}'.format(r))
print('预报和准确值的相对误差为：{:.6}%,绝对误差为:{:.6}'.format(abs(y-r)/r*100,abs(y-r)))
x1=[0];y1=[1];y2=[1]
for i in range(300):
    y1.append(g(i*0.01))
    x1.append(i*0.01)
    y2.append(no(x0,0.01+i*0.01))
plt.plot(x1,y1,c='r',label='Value')
plt.plot(x1,y2,c='k',label='prediction')
plt.xlabel('X');plt.ylabel('Y');plt.title('P-M trace')
plt.legend()
plt.show()
def f(x,y):#question second with R-K.way
    xv=mh.sin(x)+mh.cos(y)
    return xv
c1=[];c2=[]
for j in range(1000):
    x1=0;x2=10;h=0.01;y0=0;x0=0
    t=(x2-x1)/h;hi=int(t)-j
    for i in range(hi):
        n1=f(x0,y0)
        n2=f(x0+h/2,y0+n1*h/2)
        n3=f(x0+h/2,y0+n2*h/2)
        n4=f(x0+h,y0+n3*h)
        y=y0+h/6*(n1+2*n2+2*n3+n4)
        y0=y
        x0=x0+h
    c2.append(y);c1.append(hi*h)
print('问题2:(保留6位小数)\n求解X={}时,y值\nR-K求解：{:.6};步长为0.01'.format(c1[900],c2[900]))
plt.plot(c1,c2,c='k')
plt.xlabel('X');plt.ylabel('Y');plt.title('R-K trace')
plt.show()

---

结果： 

问题1:(保留6位小数)
求解X=2时,y值
预报校正法求解：0.200007;步长为0.01
准确求解微分方程：0.2
预报和准确值的相对误差为：0.00325958%,绝对误差为:6.51915e-06
问题2:(保留6位小数)
求解X=1.2307313059275027时,y值
R-K求解：1.23073;步长为0.01

---