import numpy as np
from scipy.integrate import solve_ivp
# from scipy.integrate import odeint
import matplotlib.pyplot as plt
import pandas as pd

# x0, y0: 初期条件，y0 = y'(x0)
# xs: yを求めるxの値
# dfunc: 導関数
# n: 区間数
# y: 数値解，y = y(xs)
def euler(dfunc, x0, y0, xs, n): # オイラー法（Euler method）
    h = (xs-x0)/n
    y = y0
    for i in range(n):
        x = x0+h*i
        y = y+h*dfunc(x,y)
    return y

# x0, y0: 初期条件，y0 = y'(x0)
# xs: yを求めるxの値
# dfunc: 導関数
# n: 区間数
# y: 数値解，y = y(xs)
def rk4(dfunc, x0, y0, xs, n): # 4次のルンゲ・クッタ法
    h = (xs-x0)/n
    h2 = h*0.5
    y = y0
    for i in range(n):
       x = x0+h*i
       k1 = h*dfunc(x,y)
       k2 = h*dfunc(x+h2,y+k1*0.5)
       k3 = h*dfunc(x+h2,y+k2*0.5)
       k4 = h*dfunc(x+h,y+k3)
       y = y+(k1+2.0*(k2+k3)+k4)/6.0
    return y

def dfunc7(x,y):# 微分方程式 y' = dy/dx = f(x,y)
    dydx = -4.0*(x-1.0)*y
    return dydx

def func7(x):# 微分方程式を解析的に解いて得られた関数
    y = np.exp(-2.0*(x-1.0)*(x-1.0))
    return y

eps = 1.0e-9
n, m = 500, 31
x00, y00 = 0.0, np.exp(-2.0)# 初期条件 for func7
xe = 3.0
dx = (xe-x00)/(m-1)
x0 = x00
y0a, y0b, y0c = y00, y00, y00
y_ini = np.array([y00])
print(f"{'x':>8s} {'solve_ivp':>10s} {'er':>13s} {'Euler':>10s} "\
      f"{'er2':>13s} {'rk4':>10s} {'er4':>13s}")
df = pd.DataFrame({'$x$':[],'solve_ivp':[],'er':[],'Euler':[],'er2':[],'rk4':[],'er4':[],})
for i in range(m):
    x = x00+dx*i
    y = func7(x)
    y1 = euler(dfunc7, x0, y0a, x, n)
    y4 = rk4(dfunc7, x0, y0c, x, n)
    solver = solve_ivp(dfunc7, [x0, x], y_ini, dense_output=True, method='DOP853')
    y5 = solver.sol(x)
    print(f"{x:8.3f} {y5[0]:10.6f} {y5[0]-y:13.10f} {y1:10.6f} "\
          f"{y1-y:13.10f} {y4:10.6f} {y4-y:13.10f}")
    x0 = x
    y0a, y0c = y1, y4
    y_ini = y5.copy()
    df1 = pd.DataFrame({'$x$':[x,],'solve_ivp':[y5[0],],'er':[y5[0]-y,],
                                      'Euler':[y1,],'er2':[y1-y,],
                                      'rk4':[y4,],'er4':[y4-y,]},index=[i,])
    df = pd.concat([df, df1])

       x  solve_ivp            er      Euler           er2        rk4           er4
   0.000   0.135335  0.0000000000   0.135335  0.0000000000   0.135335  0.0000000000
   0.100   0.197899  0.0000000000   0.197878 -0.0000206774   0.197899 -0.0000000000
   0.200   0.278037  0.0000000000   0.277987 -0.0000500973   0.278037 -0.0000000000
   0.300   0.375311  0.0000000000   0.375225 -0.0000864315   0.375311 -0.0000000000
   0.400   0.486752  0.0000000000   0.486627 -0.0001255903   0.486752  0.0000000000
   0.500   0.606531  0.0000000000   0.606369 -0.0001616715   0.606531  0.0000000000
   0.600   0.726149  0.0000000000   0.725961 -0.0001881401   0.726149  0.0000000000
   0.700   0.835270  0.0000000000   0.835071 -0.0001994950   0.835270  0.0000000000
   0.800   0.923116  0.0000000000   0.922923 -0.0001929179   0.923116  0.0000000000
   0.900   0.980199  0.0000000000   0.980029 -0.0001693099   0.980199  0.0000000000
   1.000   1.000000  0.0000000000   0.999867 -0.0001332711   1.000000  0.0000000000
   1.100   0.980199  0.0000000000   0.980107 -0.0000919499   0.980199  0.0000000000
   1.200   0.923116  0.0000000000   0.923063 -0.0000531152   0.923116  0.0000000000
   1.300   0.835270  0.0000000000   0.835247 -0.0000231094   0.835270  0.0000000000
   1.400   0.726149  0.0000000000   0.726144 -0.0000053686   0.726149  0.0000000000
   1.500   0.606531  0.0000000000   0.606531  0.0000000485   0.606531  0.0000000000
   1.600   0.486752  0.0000000000   0.486748 -0.0000041131   0.486752  0.0000000000
   1.700   0.375311  0.0000000000   0.375298 -0.0000135808   0.375311 -0.0000000000
   1.800   0.278037  0.0000000000   0.278013 -0.0000240024   0.278037 -0.0000000000
   1.900   0.197899  0.0000000000   0.197867 -0.0000320749   0.197899 -0.0000000000
   2.000   0.135335  0.0000000000   0.135299 -0.0000360846   0.135335 -0.0000000000
   2.100   0.088922 -0.0000000000   0.088886 -0.0000358518   0.088922  0.0000000000
   2.200   0.056135 -0.0000000000   0.056102 -0.0000322734   0.056135  0.0000000000
   2.300   0.034047 -0.0000000000   0.034021 -0.0000267288   0.034047  0.0000000000
   2.400   0.019841 -0.0000000000   0.019821 -0.0000205696   0.019841  0.0000000000
   2.500   0.011109 -0.0000000000   0.011094 -0.0000148095   0.011109  0.0000000000
   2.600   0.005976  0.0000000000   0.005966 -0.0000100241   0.005976  0.0000000000
   2.700   0.003089  0.0000000000   0.003082 -0.0000064020   0.003089  0.0000000000
   2.800   0.001534  0.0000000000   0.001530 -0.0000038686   0.001534  0.0000000000
   2.900   0.000732  0.0000000000   0.000730 -0.0000022166   0.000732  0.0000000000
   3.000   0.000335  0.0000000000   0.000334 -0.0000012064   0.000335  0.0000000000

df.style.format({"$x$":"{:.1f}", "er1":"{:.11f}", "er2":"{:.11f}", "er4":"{:.11f}"})

#solve_ivp?

plt.figure()
plt.plot(df["$x$"], df["solve_ivp"], label='solve_ivp (DOP853)')
plt.plot(df["$x$"], df["Euler"], 'o--', label='Euler')
plt.plot(df["$x$"], df["rk4"], 's--', label='RK4')
plt.xlabel("x")
plt.ylabel("y")
plt.title("Euler, RK4, solve_ivp")
plt.legend()
plt.grid()
plt.show()

	$x$	solve_ivp	er	Euler	er2	rk4	er4
0	0.0	0.135335	0.000000	0.135335	0.00000000000	0.135335	0.00000000000
1	0.1	0.197899	0.000000	0.197878	-0.00002067743	0.197899	-0.00000000000
2	0.2	0.278037	0.000000	0.277987	-0.00005009735	0.278037	-0.00000000000
3	0.3	0.375311	0.000000	0.375225	-0.00008643152	0.375311	-0.00000000000
4	0.4	0.486752	0.000000	0.486627	-0.00012559025	0.486752	0.00000000000
5	0.5	0.606531	0.000000	0.606369	-0.00016167145	0.606531	0.00000000000
6	0.6	0.726149	0.000000	0.725961	-0.00018814012	0.726149	0.00000000000
7	0.7	0.835270	0.000000	0.835071	-0.00019949496	0.835270	0.00000000000
8	0.8	0.923116	0.000000	0.922923	-0.00019291789	0.923116	0.00000000000
9	0.9	0.980199	0.000000	0.980029	-0.00016930991	0.980199	0.00000000000
10	1.0	1.000000	0.000000	0.999867	-0.00013327115	1.000000	0.00000000000
11	1.1	0.980199	0.000000	0.980107	-0.00009194987	0.980199	0.00000000000
12	1.2	0.923116	0.000000	0.923063	-0.00005311517	0.923116	0.00000000000
13	1.3	0.835270	0.000000	0.835247	-0.00002310940	0.835270	0.00000000000
14	1.4	0.726149	0.000000	0.726144	-0.00000536856	0.726149	0.00000000000
15	1.5	0.606531	0.000000	0.606531	0.00000004851	0.606531	0.00000000000
16	1.6	0.486752	0.000000	0.486748	-0.00000411307	0.486752	0.00000000000
17	1.7	0.375311	0.000000	0.375298	-0.00001358083	0.375311	-0.00000000000
18	1.8	0.278037	0.000000	0.278013	-0.00002400239	0.278037	-0.00000000000
19	1.9	0.197899	0.000000	0.197867	-0.00003207488	0.197899	-0.00000000000
20	2.0	0.135335	0.000000	0.135299	-0.00003608460	0.135335	-0.00000000000
21	2.1	0.088922	-0.000000	0.088886	-0.00003585179	0.088922	0.00000000000
22	2.2	0.056135	-0.000000	0.056102	-0.00003227338	0.056135	0.00000000000
23	2.3	0.034047	-0.000000	0.034021	-0.00002672876	0.034047	0.00000000000
24	2.4	0.019841	-0.000000	0.019821	-0.00002056956	0.019841	0.00000000000
25	2.5	0.011109	-0.000000	0.011094	-0.00001480953	0.011109	0.00000000000
26	2.6	0.005976	0.000000	0.005966	-0.00001002415	0.005976	0.00000000000
27	2.7	0.003089	0.000000	0.003082	-0.00000640203	0.003089	0.00000000000
28	2.8	0.001534	0.000000	0.001530	-0.00000386860	0.001534	0.00000000000
29	2.9	0.000732	0.000000	0.000730	-0.00000221665	0.000732	0.00000000000
30	3.0	0.000335	0.000000	0.000334	-0.00000120640	0.000335	0.00000000000

常微分方程式¶