佛山科学技术学院
实  验 报  告
课程名称                    数值分析                     
实验项目            常微分方程问题初值问题数值解法         
专业班级                姓 名         学 号            
指导教师    陈剑    成 绩           日 期          
                                                                           
一. 实验目的
1、 理解如何在计算机上实现用Euler法、改进Euler法、Runge－Kutta算法求一阶常微分方程初值问题
    的数值解。
2、 利用图形直观分析近似解和准确解之间的误差。
二、实验要求
（1） 按照题目要求完成实验内容；
（2） 写出相应的Matlab 程序；
（3） 给出实验结果(可以用表格展示实验结果)；
（4） 分析和讨论实验结果并提出可能的优化实验。
（5） 写出实验报告。
三、实验步骤
1、用Matlab编写解常微分方程初值问题的Euler法、改进Euler法和经典的四阶Runge-Kutta法。
2、给定初值问题
要求：（a）用Euler法和改进的Euler法(步长均取h=0.05)及经典的四阶Runge-Kutta法(h=0.1)求（1）的数值解，并打印的值。
(b) 用经典的四阶Runge-Kutta方法解(2)，步长分别取h=0.1, 0.05,0.025,计算并打印个点的值，与准确解比较，并列表写出在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;
matlab求导
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;
步长0.1：
x        误差
x0 = 0;
xn = 1;
y0 = 1/3;
%步长0.1：
h  = 0.1;
x  = [x0:h:xn]';
y  = R_Kutta4(x0,xn,y0,h);
步长0.05：
x        误差
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);
步长0.025：
x        误差
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和误差的关系，在不出现错误选择的情况下，步长越小，误差是越小的，而且是指数型减小。
另外通过操作可以看到三种步长下的插值和准确值函数图像：
注意第一幅图有科学计数法标注（10^10）
下面两幅图可以看出，当步长减小时到0.025时，基本就与原值一致了：
六、改进实验建议