當前位置：首頁 > 编程资源 > 编程问答 >内容正文

编程问答

常微分方程的数值解法

發布時間：2023/12/9 编程问答 27 豆豆

生活随笔收集整理的這篇文章主要介紹了常微分方程的数值解法小編覺得挺不錯的,現在分享給大家,幫大家做個參考.

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）