佛山科学技术学院实 验 报 告课程名称 数值分析 实验项目 常微分方程问题初值问题数值解法 专业班级 姓 名 学 号 指导教师 陈剑 成 绩 日 期一. 实验目的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 xy 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 个点的值，与准确解25031)(x e x y x +=-比较，并列表写出在x=0.2,0.5,0.8处，对于不同步长h 下的误差，讨论同一节点处，误差随步长的变化规律。
（c ）用Matlab 绘图函数绘制（2）的精确解和近似解的图形。
四、实验结果 %Euler.mfunction 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+1y (i) = y(i-1)+h*y_(i-1); y_(i) = f_xy(x(i),y(i)); end%Imp_Euler.mfunction 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+1y_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.mfunction 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+1y(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.mfunction y_=f_xy(x,y)%求解第五次作业第一题的点（x,y）处的导数；y_ = 1/(x^2) - y/x;%run521.mclc,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.mfunction y_=f_xy(x,y) %求第二个方程的导数： y_ = -50*y+50*(x^2)+2*x;%run522.mclc,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时，误差就很小了，再缩小一半，误差就更小了。
