微分方程数值解法课程设计许.docx

微分方程数值解法课程设计----A组

计算022班3号 许**

用复化梯形计算积分解：

算法思想：

?10

0

e?x2dx

根据复化梯形公式:取步长h?(b?a)/n，x ?a?ih(i?0,1, ,n)

i

有I(f)??b

a

N=100;a=0;b=10;

h=(b-a)/N;

f(x)dx?T

n

(f)??n

i?1

h[f(x2

i?1

)?f(x

i

)]，写出程序，如下所示：

T=zeros(0,N);forj=1:N

T(j)=0.5*h*[exp(-(a+h*(j-1)).^2)+exp(-(a+h*j).^2)]t=sum(T)

end

运行结果：t=0.8862

用Euler、改进的Euler法、梯形法、R-K法解

?y??x?y

??y（0）?1

?

（0?x?1）

y??x?1?2ex

解：

根据算法思想：

Euler法

y ?y ?h(x ?y)

n?1 n n n

改进的Euler法

?y

? n?1

?y ?h(x

n n

y)

n

?

??y

n?1

?y ?

n

h?(x

2 n

?y)?(x

n

n?1

?y )?

n?1

梯形法

y

n?1

?y ?

n

h?(x

2 n

?y)?(x

n

n?1

?y )?

n?1

（4）R-K法

y

n?1

?y ?

n

h?k

6 1

2k

2

2k

3

k?

4

?k ?f(x,y)

1?1 n n

1

?k?f(x

? 2 n

?

h,y

2 n

)2 1

1我们可以写出程序如下所示：

1

?k ?f(x

?3 n

h,y

2 n

hk)

2 2

k?? ?f(x

k

?

4 n

h,y

n

hk)

3

formatlong;N=100;

h=1/N;%步长

U=zeros(1,N+1);%eulerT=zeros(1,N+1);%真解X=zeros(1,N+1);R=zeros(4,N+1);%误差G=zeros(1,N+1);%改进euleO=zeros(1,N+1);%梯形法K=zeros(1,N+1)%R-K法c=zeros(1,4);

a=0;

forj=2:N+1;

X(j)=h*(j-1);

endU(1)=1T(1)=1;

G(1)=1;

X(1)=0;

O(1)=1;

K(1)=1;

forj=1:N

U(j+1)=U(j)+h*(U(j)+X(j));%euler法公式

a=G(j)+h*(G(j)+X(j));

G(j+1)=G(j)+h/2*(G(j)+X(j)+a+X(j+1)); %改进的euler法公式b=O(j)+h*(O(j)+X(j));

fori=1:10

O(j+1)=O(j)+h/2*(G(j)+X(j)+b+X(j+1));%梯形法公式b=O(j+1);

endc(1)=K(j)+X(j);

c(2)=K(j)+h/2*c(1)+(X(j)+h/2);

c(3)=K(j)+h/2*c(2)+(X(j)+h/2);

c(4)=K(j)+h*c(3)+(X(j)+h);

K(j+1)=K(j)+h/6*(c(1)+2*c(2)+2*c(3)+c(4));%四阶R-K法公式T(j+1)=-X(j+1)-1+2*exp(X(j+1));%真实值R(1,j+1)=U(j+1)-T(j+1);%EULER法误差R(2,j+1)=G(j+1)-T(j+1);%改进的EULER法误差R(3,j+1)=O(j+1)-T(j+1);%梯形法误差

R(4,j+1)=K(j+1)-T(j+1);%四阶R-K法误差

end

真实值T：

1 1.0101 1.0204 1.0309 1.0416 1.0525 1.0637 1.075 1.0866 1.0983 1.1103 1.1226 1.135

1.1477

1.1605

1.1737

1.187

1.2006

1.2144

1.2285

1.2428

1.2574

1.2722

1.2872

1.3025

1.3181

1.3339

1.3499

1.3663

1.3829

1.3997

1.4169

1.4343

1.4519

1.4699

1.4881

1.5067

1.5255

1.5446

1.564

1.5836

1.6036

1.6239

1.6445

1.6654

1.6