MATLAB数值分析实验五(欧拉法,荣格-库塔法解常微分方程)

佛山科学技术学院
实 验 报 告
课程名称 数值分析 实验项目 常微分方程问题初值问题数值解法 专业班级 姓 名 学 号 指导教师 陈剑 成 绩 日 期
一. 实验目的
1、理解如何在计算机上实现用Euler 法、改进Euler 法、Runge －Kutta 算法求一阶常微分方程初值问题
⎩⎨
⎧=∈='1
)(]
,[),,()(y a y b a x y x f x y 的数值解。

2、利用图形直观分析近似解和准确解之间的误差。

二、实验要求
（1） 按照题目要求完成实验内容； （2） 写出相应的Matlab 程序；
（3） 给出实验结果(可以用表格展示实验结果)； （4） 分析和讨论实验结果并提出可能的优化实验。

（5） 写出实验报告。

三、实验步骤
1、用Matlab 编写解常微分方程初值问题的Euler 法、改进Euler 法和经典的四阶Runge-Kutta 法。

2、给定初值问题
⎪⎩⎪⎨⎧
=≤≤-=;
1)1(,
21,1')1(2y x x
y x y
⎪
⎩
⎪⎨⎧=≤≤++-=31)0(10,25050')2(2y x x x y y 要求：（a ）用Euler 法和改进的Euler 法(步长均取h=0.05)及经典的四阶Runge-Kutta 法(h=0.1)求（1）的数值解，并打印)10,....2,1,0(1.01=+=i i x 的值。

(b) 用经典的四阶Runge-Kutta 方法解(2)，步长分别取h=0.1, 0.05,0.025,计算并打印
)10,....2,1,0(1.0==i i x 个点的值，与准确解2503
1
)(x e x y x +=-比较，并列表写出在
x=0.2,0.5,0.8处，对于不同步长h 下的误差，讨论同一节点处，误差随步长的变化规律。

（c ）用Matlab 绘图函数绘制（2）的精确解和近似解的图形。

四、实验结果 %Euler.m
function y = Euler(x0,xn,y0,h) %Euler 法解方程f_xy ； %x0,y0为初始条件； %x0,xn 为求值区间； %h 为步长； %求区间个数： n = (xn-x0)/h;
%矩阵x 存储n+1个节点： x = [x0:h:xn]';
%矩阵y 存储节点处的值： y = [y0;zeros(n,1)];
%矩阵y_存储节点处导数值： y_(1)= f_xy(x0,y0); y_ = [y_(1);zeros(n,1)];
%进行迭代（欧拉法迭代；求导数）： for i = 2:n+1
y (i) = y(i-1)+h*y_(i-1); y_(i) = f_xy(x(i),y(i)); end
%Imp_Euler.m
function y = Imp_Euler(x0,xn,y0,h)
%改进的Euler法解方程f_xy；
%x0,y0为初始条件；
%x0,xn为求值区间；
%h为步长；
%求区间个数：
n = (xn-x0)/h;
%矩阵x存储n+1个节点：
x = [x0:h:xn]';
%矩阵y存储节点处的值：
y = [y0;zeros(n,1)];
%矩阵y_存储节点处导数值：
y_(1)= f_xy(x0,y0);
y_ = [y_(1);zeros(n,1)];
%使用改进Euler法求值（欧拉法求近似；近似点导数；梯形校正；求导）：for i = 2:n+1
y_l = y(i-1) + h*y_(i-1);
y_l = f_xy(x(i),y_l);
y(i) = y(i-1) + (h/2)*(y_(i-1)+y_l);
y_(i)= f_xy(x(i),y(i));
end
%R_Kutta4.m
function y = R_Kutta4(x0,xn,y0,h)
%Runger_Kutta法解方程f_xy；
%x0,y0为初始条件；
%x0,xn为求值区间；
%h为步长；
%求区间个数：
n = (xn-x0)/h;
%矩阵x存储n+1个节点：
x = [x0:h:xn]';
%矩阵y存储节点处的值：
y = [y0;zeros(n,1)];
%矩阵k1,k2,k3,k4存储各节点（中点）数值：
k1(1)= f_xy(x0,y0);
k1 = [k1(1);zeros(n,1)];
k2(1)= f_xy(x0+h/2,y0+h*k1(1)/2);
k2 = [k2(1);zeros(n,1)];
k3(1)= f_xy(x0+h/2,y0+h*k2(1)/2);
k3 = [k3(1);zeros(n,1)];
k4(1)= f_xy(x0+h,y0+h*k3(1));
k4 = [k4(1);zeros(n,1)];
for i= 2:n+1
y(i) = y(i-1)+(h/6)*(k1(i-1)+2*k2(i-1)+2*k3(i-1)+k4(i-1));
k1(i)= f_xy(x(i),y(i));
k2(i)= f_xy(x(i)+h/2,y(i)+h*k1(i)/2);
k3(i)= f_xy(x(i)+h/2,y(i)+h*k2(i)/2);
k4(i)= f_xy(x(i)+h,y(i)+h*k3(i));
end
（a）：
%f_xy.m
function y_=f_xy(x,y)
%求解第五次作业第一题的点（x,y）处的导数；
y_ = 1/(x^2) - y/x;
%run521.m
clc,clear;
x0 = 1;
xn = 2;
h = 0.05;
y0 = 1;
%便于显示出x，与y对应：x = [x0:h:xn]';
y = Euler(x0,xn,y0,h);
YE =[x,y];
y = Imp_Euler(x0,xn,y0,h); YIE= [x,y];
h = 0.1;
x = [x0:h:xn]';
y = R_Kutta4(x0,xn,y0,h); YRK= [x,y];
(b)： %f_xy.m
function y_=f_xy(x,y) %求第二个方程的导数： y_ = -50*y+50*(x^2)+2*x;
%run522.m
clc,clear; x0 = 0; xn = 1; y0 = 1/3; %步长0.1： h = 0.1; x = [x0:h:xn]';
y = R_Kutta4(x0,xn,y0,h); y_r= f_Real(x); Y1 = [x,y,y_r];
%步长0.05： h = 0.05; x = [x0:h:xn]';
y = R_Kutta4(x0,xn,y0,h); y_r= f_Real(x); Y2 = [x,y,y_r]; %步长0.025： h = 0.025; x = [x0:h:xn]';
y = R_Kutta4(x0,xn,y0,h); y_r= f_Real(x); Y3 = [x,y,y_r];
五、讨论分析
(a)
从结果可以看出使用RK 方法，步长较大但是结果也更加精确； (b)
分析求值结果的误差，可以发现当步长取0.1时，误差是超级大的（10^8数量级），但是当步长缩小一半取0.05时，误差就很小了，再缩小一半，误差就更小了。

根据书本的介绍，步长缩小一半，精度提高十多倍，但是从测试情况来看，步长有个临界点（可能仅仅针对这种凹形，指数形函数），太大时，会忽略掉函数的某些属性，譬如取0.1时跳过了凹形区域，从而丧失了精确度（错误结果），所以步长选择时要比较谨慎，以免出现这种情况，至于h 和误差的关系，在不出现
步长0.1： x 误差
步长0.05： x 误差
步长0.025： x 误差
错误选择的情况下，步长越小，误差是越小的，而且是指数型减小。

另外通过操作可以看到三种步长下的插值和准确值函数图像：
注意第一幅图有科学计数法标注（10^10）
下面两幅图可以看出，当步长减小时到0.025时，基本就与原值一致了：
六、改进实验建议
实验的题目表述有比较大的问题，有重复操作部分，有些表达也不是很清楚，所以希望出题时能够认真考虑下实现过程，是否有矛盾或是冗余的，不然很浪费时间又不能锻炼操作能力。

总的来说还是有收获的，尤其是RK步长取0.1时都觉得是否程序出错。