实验九 欧拉方法
信息与计算科学金融 崔振威 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(,22'ttetyyytyt
解:
执行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 = 0 0.0500 0.1000 0.1500 0.2000 0.2500 0.3000 0.3500
0.4000 0.4500 0.5000 0.5500 0.6000 0.6500 0.7000 0.7500
0.8000 0.8500 0.9000 0.9500 1.0000 1.0500 1.1000 1.1500
1.2000 1.2500 1.3000 1.3500 1.4000 1.4500 1.5000 1.5500
1.6000 1.6500 1.7000 1.7500 1.8000 1.8500 1.9000 1.9500
2.0000
y1 =
1.0000 0.9500 0.9026 0.8580 0.8162 0.7774 0.7417 0.7091
0.6798 0.6538 0.6312 0.6121 0.5967 0.5848 0.5767 0.5724
0.5719 0.5753 0.5826 0.5940 0.6094 0.6290 0.6526 0.6805
0.7126 0.7490 0.7897 0.8347 0.8841 0.9379 0.9961 1.0588
1.1260 1.1977 1.2739 1.3547 1.4401 1.5301 1.6247 1.7240
1.8279
在Matlab 命令窗口输入：
>> y3=-exp(-x1)+x1.^2-2*x1+2
求得解析解：
y3 =
1.0000 0.9513 0.9052 0.8618 0.8213 0.7837 0.7492 0.7178
0.6897 0.6649 0.6435 0.6256 0.6112 0.6005 0.5934 0.5901
0.5907 0.5951 0.6034 0.6158 0.6321 0.6526 0.6771 0.7059
0.7388 0.7760 0.8175 0.8633 0.9134 0.9679 1.0269 1.0903
1.1581 1.2305 1.3073 1.3887 1.4747 1.5653 1.6604 1.7602
1.8647
输入：
>> plot(x1,y1,'o');hold on
>> plot(x1,y3,'*');hold on
可得近似值与解析解的图像比较：
