资源描述:
《仿真建模作业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