常微分方程的数值解法

1 实验环境

Matlab、pycharm环境

2 实验目的

掌握求解常微分方程的欧拉法，Runge-Kutta方法（二阶的预估校正法，经典的四阶R-K法）

3实验原理

4实验内容

1）实验方案

方案一：用欧拉法，预估校正法，经典的四阶龙格库塔方法求解下列ODE问题：

例题：在区间【0，1】上以h=0.1用欧拉法，预估校正法,经典的四阶龙格库塔法求解微分方程 dy/dx=-y+x+1，初值y（0）=1；其精确解为y=x+exp(-x),且将计算结果与精确解进行比较，对三个算法的收敛性的进行分析比较。

方案二：用欧拉法，预估校正法, 经典的四阶龙格库塔方法求解初值问题

2）实验步骤

用matlab编写欧拉法和预估校正法程序

用python编写经典四阶龙格库塔程序

数据处理

方案一（h=0.05)
n	xn	yn
n	xn	欧拉法	预估校正法	经典四阶龙格库塔法	准确解
0	0	1.0000	1.0000	1.00000000	1.00000000
1	0.05	1.0000	1.0012	1.00122943	1.00122942
2	0.1	1.0025	1.0049	1.00483742	1.00483742
3	0.15	1.0074	1.0108	1.01070798	1.01070798
4	0.2	1.0145	1.0188	1.01873076	1.01873075
5	0.25	1.0238	1.0289	1.02880079	1.02880078
6	0.3	1.0351	1.0409	1.04081823	1.04081822
7	0.35	1.0483	1.0548	1.05468810	1.05468809
8	0.4	1.0634	1.0704	1.07032006	1.07032005
9	0.45	1.0802	1.0878	1.08762817	1.08762815
10	0.5	1.0987	1.1067	1.10653068	1.10653066
11	0.55	1.1188	1.1271	1.12694983	1.12694981
12	0.6	1.1404	1.1490	1.14881165	1.14881164
13	0.65	1.1633	1.1722	1.17204580	1.17204578
14	0.7	1.1877	1.1967	1.19658532	1.19658530
15	0.75	1.2133	1.2225	1.22236657	1.22236655
16	0.8	1.2401	1.2495	1.24932898	1.24932896
17	0.85	1.2681	1.2776	1.27741495	1.27741493
18	0.9	1.2972	1.3067	1.30656968	1.30656966
19	0.95	1.3274	1.3369	1.33674104	1.33674102
20	1	1.3585	1.3680	1.36787946	1.36787944

方案二(h=0.05)
n	xn	yn
n	xn	欧拉法	预估校正法	经典四阶龙格库塔法	准确解
0	0	1	1	1.00000000	1.00000000
1	0.05	0.95	0.9524	0.95241847	0.95241846
2	0.1	0.9049	0.9095	0.90936161	0.90936161
3	0.15	0.8642	0.8708	0.87039095	0.87039094
4	0.2	0.8274	0.8358	0.83510538	0.83510537
5	0.25	0.7942	0.8042	0.80313832	0.80313831
6	0.3	0.7642	0.7757	0.77415506	0.77415504
7	0.35	0.7371	0.7499	0.74785026	0.74785024
8	0.4	0.7126	0.7265	0.72394567	0.72394565
9	0.45	0.6904	0.7053	0.70218802	0.70218800
10	0.5	0.6702	0.686	0.68234702	0.68234699
11	0.55	0.6519	0.6685	0.66421349	0.66421347
12	0.6	0.6351	0.6525	0.64759775	0.64759773
13	0.65	0.6199	0.6378	0.63232797	0.63232795
14	0.7	0.6058	0.6244	0.61824873	0.61824870
15	0.75	0.5929	0.612	0.60521967	0.60521965
16	0.8	0.581	0.6004	0.59311426	0.59311423
17	0.85	0.5699	0.5897	0.58181860	0.58181858
18	0.9	0.5596	0.5797	0.57123039	0.57123037
19	0.95	0.5499	0.5703	0.56125793	0.56125791
20	1	0.5408	0.5614	0.55181918	0.55181916

方案二(h=0.05)
n	xn	yn
n	xn	欧拉法	预估校正法	经典四阶龙格库塔法	准确解
0	0	1	1	1.00000000	1.00000000
1	0.1	0.9	0.9095	0.90936175	0.90936161
2	0.2	0.819	0.8363	0.83510562	0.83510537
3	0.3	0.7535	0.777	0.77415537	0.77415504
4	0.4	0.7004	0.7288	0.72394602	0.72394565
5	0.5	0.6572	0.6895	0.68234739	0.68234699
6	0.6	0.6218	0.6572	0.64759814	0.64759773
7	0.7	0.5925	0.6303	0.61824911	0.61824870
8	0.8	0.568	0.6076	0.59311463	0.59311423
9	0.9	0.5472	0.5881	0.57123076	0.57123037
10	1	0.5291	0.571	0.55181953	0.55181916

图像分析

方案一欧拉与预估（20）

方案一D-K（20）

方案二欧拉与预估（10）

方案二D-K（10）

方案二欧拉与预估（20）

方案二D-K（20）

5实验结论

根据实验数据和实验图像分析可知，收敛性欧拉法<预估校正法<经典四阶龙格库塔法。

附录1：源程序

方案一程序：用欧拉法，预估校正法（matlab）

function [y0,y1]=Euler(a,b,h,n)
% 精确解
t=0:0.01:1;
y=t+exp(-t);    %精确解
plot(t,y,'b'); hold on;
% 初始化
a=0.0;b=1.0;n=20;h=(b-a)/n; t0=a:h:b;
%  Euler’s method
y0(1)=1;  %初值
for k=1:20
    y0(k+1)=y0(k)+h*(-y0(k)+t0(k)+1);     %y0输出欧拉法解
end
plot(t0,y0,'ro'); hold on;
%  predictor-corrector method
y1(1)=1.0;
for i=1:20
    y1(i+1)=y1(i)+h*(-y1(i)+t0(i)+1);
    y1(i+1)=y1(i)+h*(-y1(i)+t0(i)-y1(i+1)+t0(i+1)+2)/2;   %y1输出预估校正法解
end
plot(t0,y1,'bs')

方案一程序：D-K（python）

import numpy as np

import math

import matplotlib.pyplot as plt

def fun(i):

return i+math.exp(-i)

def DK(a,b,n):

'''绘制精确解'''

h=(b-a)/n

#精确解

x=np.arange(0,1+0.01,0.01)

y=[]

for i in x:

y.append(fun(i))

'''绘制D-K散点图'''

x0=np.arange(0,1+h,h)

y0=[1]

for i in range(n):

k1 = h * (-y0[i] + x0[i] + 1)

k2 = h * (-y0[i]-0.5*k1 + x0[i]+0.5*h + 1)

k3 = h * (-y0[i]-0.5*k2 + x0[i]+0.5*h + 1)

k4 = h * (-y0[i]-k3 + x0[i]+h + 1)

y0.append(y0[i]+(k1+2*k2+2*k3+k4)/6)

for i in y0:

print('{:.4f}'.format(i),end=';')

plt.title('i+math.exp(-i)')

plt.plot(x, y, color='skyblue', label='true')

plt.scatter(x0, y0, c='c', label='D-K')

plt.legend()

plt.xlabel('x')

plt.ylabel('y')

plt.show()

DK(0,1,20)

方案二程序：用欧拉法，预估校正法（matlab）

function [y0,y1]=Euler(a,b,h,n)
% 精确解
t=0:0.01:1;
y=0.5.*(t.^2+2).*exp(-t);    %精确解
plot(t,y,'b'); hold on;
% 初始化
a=0.0;b=1.0;n=10;h=(b-a)/n; t0=a:h:b;
%  Euler’s method
y0(1)=1;  %初值
for k=1:10
    y0(k+1)=y0(k)+h*(t0(k)*exp(-t0(k))-y0(k));     %y0输出欧拉法解
end
plot(t0,y0,'ro'); hold on;
%  predictor-corrector method
y1(1)=1.0;
for i=1:10
    y1(i+1)=y1(i)+h*(t0(i)*exp(-t0(i))-y0(i));
    y1(i+1)=y1(i)+h*(t0(i)*exp(-t0(i))-y0(i)+t0(i+1)*exp(-t0(i+1))-y0(i+1))/2;   %y1输出预估校正法解
en
plot(t0,y1,'bs')

方案二程序：用D-K（python）