常微分方程初值问题数值解法

《常微分方程初值问题数值解法》由会员上传分享，免费在线阅读，更多相关内容在教育资源-天天文库。

1、贵州师范大学数学与计算机科学学院学生实验报告课程名称：数值分析班级实验日期：学号：姓名指导教师：实验成绩：一、实验名称实验六:常微分方程初值问题数值解法二、实验目的及要求1.让学生掌握用Euler法,Runge-Kutta法求解常微分方程初值问题.2.培养Matlab编程与上机调试能力.三、实验环境每人一台计算机，要求安装WindowsXP操作系统，Microsoftoffice2003、MATLAB6.5(或7.0).四、实验内容1.取步长h=0.1,0.05,0.01,,用Euler法及经典4阶Runge-Kutta法求解初值问题要求:1)画出准确解(准确解)的曲线,近似解折线;

2、2)把节点0.1和0.5上的精确解与近似解比较,观察误差变化情况.2.用Euler法,隐式Euler法和经典4阶R-K法取不同步长解初值问题并画出曲线观察稳定性.注：题1必须写实验报告五、算法描述及实验步骤Euler法：输入输出Euler解y步1步2对执行步3输出经典4阶R-K法：输入输出4阶R-K解y步1步2对执行，，，步3输出六、调试过程及实验结果>>a=0;b=1;y0=1;h=0.1;x0=0;>>RK('fun',a,b,y0,h)'ans=0.82870.71030.63880.60930.61790.66120.73660.84190.97531.1353>>Euler

3、('fun',a,b,x0,y0,h)ans=1.00000.80000.66200.57760.54010.54410.58530.66020.76620.90091.0627>>E('fun',a,b,x0,y0,h)err=0.02870.0738ans=1.00000.80000.66200.57760.54010.54410.58530.66020.76620.90091.0627七、总结1、在调试是常出现错误，主要是一些细节没注意，如一个“：”打成“；”等。2、对于第2）小问在画图时是在老师的指导下才有了清晰的思路，下来才晓得怎么编程。3、还存在程序中“h”我只会赋值，我

4、尝试用循环语句并调用一些源程序，但最终还是没调试出来，但我还是收获很多，也让我回忆了很多知识及其巩固了一些知识点。八、附录（源程序清单）四阶Runge-Kutta法functionY=RK(f,a,b,y0,h)m=(b-a)/h;Y=zeros(m,1);x=a;forn=1:mK1=feval(f,x,y0);x=x+0.5*h;y1=y0+0.5*h*K1;K2=feval(f,x,y1);y2=y0+0.5*h*K2;K3=feval(f,x,y2);x=x+0.5*h;y3=y0+h*K3;K4=feval(f,x,y3);y0=y0+h*(K1+2*K2+2*K3+K4)

5、/6;Y(n)=y0;end函数定义functionz=fun(x,y)z=-2*y+2*x^2+2*x;欧拉法及准确解图像Functiony=Euler(f,a,b,x0,y0,h)m=(b-a)/h;y=zeros(1,m+1);x(1)=x0;fork=1:mx(k+1)=a+k*h;endy=zeros(1,m+1);y(1)=y0;fork=1:my(k+1)=y(k)+h*feval(f,x(k),y(k));endplot(x,y,'b-')holdont=a:0.01:b;y0=exp(-2*t)+t.^2;plot(t,y0,'r-')text(t(20),y0(2

6、0),'leftarrowy=exp(-2*t)+t^2')title('׼ȷ½âÇúÏßÓë½üËƽâÕÛÏß')text(x(8),y(8),'leftarrowy=Euler½üËƽâÕÛÏß')0.1与0.5处误差及图像functiony=E(f,a,b,x0,y0,h)m=(b-a)/h;y=zeros(1,m+1);x(1)=x0;m1=(0.1-a)/h;m2=(0.5-a)/h;fork=1:mx(k+1)=a+k*h;endy=zeros(1,m+1);y(1)=y0;fork=1:my(k+1)=y(k)+h*feval(f,x(k),y(k));en

7、dy0=y;X=[x(m1+1),x(m2+1)];Y=[y0(m1+1),y0(m2+1)];t1=0.1;t2=0.5;T=[t1,t2];Y1=exp(-2*t1)+t1.^2;Y2=exp(-2*t2)+t2.^2;T1=[Y1,Y2];err=T1-Yholdonplot(x,y,'b-')holdont=a:0.01:b;y0=exp(-2*t)+t.^2;plot(t,y0,'r-')text(t(20),y0(20),'leftarrow