Runge-Kutta法

#P240  例3
import numpy as np

accuracy=6

def Fx(x):
    return x-1+np.e**(-x)

def fx(x,y):
    return x-y

x=0
y=0
h=0.2
for i in range(0,5):
    k1=fx(x,y)
    k2=fx(x+h/2,y+h/2*k1)
    k3=fx(x+h/2,y+h/2*k2)
    k4=fx(x+h,y+h*k3)
    y1=y+h*(k1+2*k2+2*k3+k4)/6
    yn=round(Fx(x),accuracy)
    err_y=round(yn-y,accuracy)
    print(x,y,yn,err_y)    
    y=round(y1,accuracy)
    x=round(x+h,accuracy)

#P264  4
import numpy as np

accuracy=5

def Fx(x):
    return x*(3+x**2)/3/(1+x**2)

def fx(x,y):
    return 1-2*x*y/(1+x**2)

x=0
y=0
h=0.5
for i in range(0,5):
    k1=fx(x,y)
    k2=fx(x+h/2,y+h/2*k1)
    k3=fx(x+h/2,y+h/2*k2)
    k4=fx(x+h,y+h*k3)
    y1=y+h*(k1+2*k2+2*k3+k4)/6
    yn=round(Fx(x),accuracy)
    err_y=round(yn-y,accuracy)
    #print(k1,k2,k3,k4,sep='\t')
    print(x,y,yn,err_y)    
    y=round(y1,accuracy)
    x=round(x+h,accuracy)

#P264  7
import numpy as np

accuracy=4

def Fx(x):
    return np.e**(-x)

def fx(x,y):
    return -y

x=0
y=1
h=0.1
for i in range(0,11):
    k1=fx(x,y)
    k2=fx(x+h/2,y+h/2*k1)
    k3=fx(x+h/2,y+h/2*k2)
    k4=fx(x+h,y+h*k3)
    y1=y+h*(k1+2*k2+2*k3+k4)/6
    yn=round(Fx(x),accuracy)
    err_y=round(yn-y,accuracy)
    #print(k1,k2,k3,k4,sep='\t')
    print(i,x,y,yn,err_y)    
    y=round(y1,accuracy)
    x=round(x+h,accuracy)

#P264  9
import numpy as np

accuracy=4

def Fx(x):
    return 2/(x**2+1)

def fx(x,y):
    return -x*y**2

x=0
y=2
h=0.2
for i in range(0,6):
    k1=fx(x,y)
    k2=fx(x+h/2,y+h/2*k1)
    k3=fx(x+h/2,y+h/2*k2)
    k4=fx(x+h,y+h*k3)
    y1=y+h*(k1+2*k2+2*k3+k4)/6
    yn=round(Fx(x),accuracy)
    err_y=round(yn-y,accuracy)
    #print(k1,k2,k3,k4,sep='\t')
    print(x,y,yn,err_y)    
    y=round(y1,accuracy)
    x=round(x+h,accuracy)