资源描述：
《仿真建模作业5》由会员上传分享，免费在线阅读，更多相关内容在学术论文-天天文库。
1、四川理工学院系统建模与仿真课程报告学生：裴俊丞学号：09021010215专业：自动化班级：自动化092任课教师：谭飞四川理工学院自动化与电子信息学院二O一一年十二月13目录第一章………………………………………………………31.1题目………………………………………………………31.2程序………………………………………………………31.3小结……………………………………………………8第二章………………………………………………………82.1题目………………………………………………………82.2解答………………………………………………………8第三章……………………………………………
2、…………123.1题目……………………………………………………123.2解答……………………………………………………12第四章………………………………………………………134.1题目……………………………………………………134.2解答……………………………………………………1313第一章1.1题目编写微分方程dy/dx＝xy,当x＝0时y＝1+5/10+9/100=1.59,x属于0～3之间，编写积分程序，包括欧拉数值积分程序，预报校正数字积分程序、4阶龙格库塔积分程序，它们的积分步长分别取0.01，0.1,0.5，绘制积分结果曲线,比较在同一步长下不同算法的误差和同一算法
3、在不同步长下的误差，得出结论说明（绿色线为欧拉法曲线，红色为预报校正法曲线，蓝色为4阶龙格库塔法曲线）1.2程序欧拉数字积分程序function[x,y]=naler(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1;y(n+1)=y(n)+h*feval(dyfun,x(n),y(n));endx=x';y=y';预报校正数字积分程序function[x,y]=zseuler(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:lengt
4、h(x)-1k1=feval(dyfun,x(n),y(n));y(n+1)=y(n)+h*k1;k2=feval(dyfun,x(n+1),y(n+1));y(n+1)=y(n)+h*(k1+k2)/2;endx=x';y=y';4阶龙格库塔积分程序function[x,y]=do(dyfun,xspan,y0,h)x=xspan(1):h:xspan(2);y(1)=y0;forn=1:length(x)-1k1=feval(dyfun,x(n),y(n));k2=feval(dyfun,x(n)+h/2,y(n)+h/2*k1);k3=feval(dyfun,x(n)+
5、h/2,y(n)+h/2*k2);k4=feval(dyfun,x(n)+h,y(n)+h*k3);13y(n+1)=y(n)+h*(k1+2*k2+2*k3+k4)/6;endx=x';y=y'步长为0.01,不同算法的积分clear;dyfun=inline('x*y');[x,y]=naler(dyfun,[0,3],1.59,0.01);figure(1)holdonplot(x,y,'g')[x,y]=zseuler(dyfun,[0,3],1.59,0.01);plot(x,y,'r')[x,y]=do(dyfun,[0,3],1.59,0.01);plot(x,
6、y)holdoffgridon图1-1步长为0.1,不同算法的积分clear;dyfun=inline('x*y');[x,y]=naler(dyfun,[0,3],1.59,0.1);figure(1)holdonplot(x,y,'g')[x,y]=zseuler(dyfun,[0,3],1.59,0.1);plot(x,y,'r')[x,y]=do(dyfun,[0,3],1.59,0.1);plot(x,y)holdoffgridon13图1-2步长为0.5,不同算法的积分clear;dyfun=inline('x*y');[x,y]=naler(dyfun,[0,3
7、],1.59,0.5);figure(2)holdonplot(x,y,'g')[x,y]=zseuler(dyfun,[0,3],1.59,0.5);plot(x,y,'r')[x,y]=do(dyfun,[0,3],1.59,0.5);plot(x,y)holdoffgridon图1-3步长分别为0.01，0.1,0.5时的欧拉法曲线13clear;dyfun=inline('x*y');[x,y]=naler(dyfun,[0,3],1.59,0.01);figure(3)holdonplot(x
