常微分方程的数值解法

方案一（h=0.05)

n

xn

yn

欧拉法

预估校正法

经典四阶龙格库塔法

准确解

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

欧拉法

预估校正法

经典四阶龙格库塔法

准确解

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

欧拉法

预估校正法

经典四阶龙格库塔法

准确解

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）

方案一程序：用欧拉法，预估校正法（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）

import numpy as np

import math

import matplotlib.pyplot as plt

def fun(i):

    return 0.5*(i**2+2)*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 * (x0[i]*math.exp(-x0[i])-y0[i])

        k2 = h * ((x0[i]+0.5*h)*math.exp(-(x0[i]+0.5*h))-y0[i]-0.5*k1)

        k3 = h * ((x0[i]+0.5*h)*math.exp(-(x0[i]+0.5*h))-y0[i]-0.5*k2)

        k4 = h * ((x0[i]+h)*math.exp(-(x0[i]+h))-y0[i]-k3)

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

    for i in y0:

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

    plt.title('0.5*(i**2+2)*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)