四阶龙格库塔和向前欧拉方法和隐式欧拉法以及改进欧拉法教材

实验九欧拉方法

信息与计算科学金融崔振威201002034031

一、实验目的：

1、掌握欧拉算法设计及程序实现

二、实验内容：

1、p364-9.2.4、p386-9.5.6

三、实验要求：

主程序：

欧拉方法（前项）：

function [x,y]=DEEuler(f,a,b,y0,n);

%f:一阶常微分方程的一般表达式的右端函数

%a:自变量的取值下限

%b:自变量的取值上限

%y0：函数的初值

%n:积分的步数

if nargin<5,n=50;

end

h=(b-a)/n;

x(1)=a;y(1)=y0;

for i=1:n

x(i+1)=x(i)+h;

y(i+1)=y(i)+h*feval(f,x(i),y(i));

end

format short

欧拉方法（后项）：

function [x,y]=BAEuler(f,a,b,y0,n);

%f:一阶常微分方程的一般表达式的右端函数

%a:自变量的取值下限

%b:自变量的取值上限

%y0：函数的初值

%n:积分的步数

if nargin<5,n=50;

end

h=(b-a)/n;

x(1)=a;y(1)=y0;

for i=1:n

x(i+1)=x(i)+h;

y1=y(i)+h*feval(f,x(i),y(i));

y(i+1)=y(i)+h*feval(f,x(i+1),y1);

end

format short

梯形算法：

function [I,step,h2] = CombineTraprl(f,a,b,eps)

%f 被积函数

%a,b 积分上下限

%eps 精度

%I 积分结果

%step 积分的子区间数

if(nargin ==3)

eps=1.0e-4;

end

n=1;

h=(b-a)/2;

I1=0;

I2=(subs(sym(f),findsym(sym(f)),a)+subs(sym(f),findsym(sym(f)),b))/h;

while abs(I2-I1)>eps

n=n+1;

h=(b-a)/n;

I1=I2;

I2=0;

for i=0:n-1

x=a+h*i;

x1=x+h;

I2=I2+(h/2)*(subs(sym(f),findsym(sym(f)),x)+subs(sym(f),findsym(sym(f)),x1));

end

end

I=I2;

step=n;

h2=(b-a)/n;

改进欧拉方法：

function [x,y]=MoEuler(f,a,b,y0,n);

%f:一阶常微分方程的一般表达式的右端函数

%a:自变量的取值下限

%b:自变量的取值上限

%y0：函数的初值

%n:积分的步数

if nargin<5,n=50;

end

h=(b-a)/n;

x(1)=a;y(1)=y0; for i=1:n

x(i+1)=x(i)+h;

y1=y(i)+h*feval(f,x(i),y(i)); y2=y(i)+h*feval(f,x(i+1),y1); y(i+1)=(y1+y2)/2; end

format short

四阶龙格-库塔法：

function y = DELGKT4_lungkuta(f, h,a,b,y0,varvec) %f:一阶常微分方程的一般表达式的右端函数 %h:积分步长

%a ：自变量的取值下限 %b:自变量的取值上限

%varvec ：常微分方程的变量组 format long; N = (b-a)/h;

y = zeros(N+1,1); y(1) = y0; x = a:h:b;

var = findsym(f); for i=2:N+1

K1 = Funval(f,varvec,[x(i-1) y(i-1)]);%Funval 为程序所需要的函数 K2 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K1*h/2]); K3 = Funval(f,varvec,[x(i-1)+h/2 y(i-1)+K2*h/2]); K4 = Funval(f,varvec,[x(i-1)+h y(i-1)+h*K3]); y(i) = y(i-1)+h*(K1+2*K2+2*K3+K4)/6; end

format short;

p364-9.2.4

欧拉方法（前项）：

1、22)(,1)0(,2

2

'

+-+-==-=-t t e t y y y t y t

解:

执行20步时：

编写函数文件doty.m 如下： function f=doty(x,y); f=x^2-y;

在Matlab 命令窗口输入：

>> [x1,y1]=DEEuler('doty',0,2,1,20)

可得到结果：

x1 =

0 0.1000 0.2000 0.3000 0.4000 0.5000 0.6000 0.7000

0.8000 0.9000 1.0000 1.1000 1.2000 1.3000 1.4000 1.5000

1.6000 1.7000 1.8000 1.9000

2.0000

y1 =

1.0000 0.9000 0.8110 0.7339 0.6695 0.6186 0.5817 0.5595

0.5526 0.5613 0.5862 0.6276 0.6858 0.7612 0.8541 0.9647

1.0932 1.2399 1.4049 1.5884 1.7906

在Matlab 命令窗口输入：

>> y3=-exp(-x1)+x1.^2-2*x1+2

求得解析解：

y3 =

1.0000 0.9052 0.8213 0.7492 0.6897 0.6435 0.6112 0.5934

0.5907 0.6034 0.6321 0.6771 0.7388 0.8175 0.9134 1.0269

1.1581 1.3073 1.4747 1.6604 1.8647

输入：

>> plot(x1,y1,'o');hold on

>> plot(x1,y3,'*');hold on

可得近似值与解析解的图像比较：

执行40步时：

在Matlab 命令窗口输入：

>> [x1,y1]=DEEuler('doty',0,2,1,40)

可得到结果：

x1 =