




已阅读5页,还剩6页未读, 继续免费阅读
版权说明:本文档由用户提供并上传,收益归属内容提供方,若内容存在侵权,请进行举报或认领
文档简介
常微分方程数值解实验报告 学院:数学与信息科学 专业:信息与计算科学 姓名:郑思义 学号:201216524 课程:常微分方程数值解实验一:常微分方程的数值解法1、 分别用Euler法、改进的Euler法(预报校正格式)和SK法求解初值问题。(h=0.1)并与真解作比较。1.1实验代码:%欧拉法function x,y=naeuler(dyfun,xspan,y0,h)%dyfun是常微分方程,xspan是x的取值范围,y0是初值,h是步长x=xspan(1):h:xspan(2);y(1)=y0; for n=1:length(x)-1 y(n+1)=y(n)+h*feval(dyfun,x(n),y(n);end %改进的欧拉法function x,m,y=naeuler2(dyfun,xspan,y0,h)%dyfun是常微分方程,xspan是x的取值范围,y0是初值,h是步长。%返回值x为x取值,m为预报解,y为校正解x=xspan(1):h:xspan(2);y(1)=y0; m=zeros(length(x)-1,1);for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n);y(n+1)=y(n)+h*k1; m(n)=y(n+1);k2=feval(dyfun,x(n+1),y(n+1);y(n+1)=y(n)+h*(k1+k2)/2;end %四阶SK法function x,y=rk(dyfun,xspan,y0,h)%dyfun是常微分方程,xspan是x的取值范围,y0是初值,h是步长。x=xspan(1):h:xspan(2);y(1)=y0; for n=1:length(x)-1 k1=feval(dyfun,x(n),y(n); k2=feval(dyfun,x(n)+h/2,y(n)+(h*k1)/2); k3=feval(dyfun,x(n)+h/2,y(n)+(h*k2)/2); k4=feval(dyfun,x(n)+h,y(n)+h*k3); y(n+1)=y(n)+(h/6)*(k1+2*k2+2*k3+k4); end%主程序x=0:0.1:1;y=exp(-x)+x;dyfun=inline(-y+x+1); x1,y1=naeuler(dyfun,0,1,1,0.1);x2,m,y2=naeuler2(dyfun,0,1,1,0.1);x3,y3=rk(dyfun,0,1,1,0.1);plot(x,y,r,x1,y1,+,x2,y2,*,x3,y3,o);xlabel(x);ylabel(y);legend(y为真解,y1为欧拉解,y2为改进欧拉解,y3为SK解,Location,NorthWest);1.2实验结果:x真解y欧拉解y1预报值m校正值y2SK解y30.0 1.00001.00001.00001.00000.1 1.00481.00001.00001.00501.00480.2 1.01871.01001.01451.01901.01870.3 1.04081.02901.03711.04121.04080.4 1.07031.05611.06711.07081.07030.5 1.10651.09051.10371.10711.10650.6 1.14881.13141.14641.14941.14880.7 1.19661.17831.19451.19721.19660.8 1.24931.23051.24751.25001.24930.9 1.30661.28741.30501.30721.30661.0 1.36791.34871.36651.36851.36792、 选取一种理论上收敛但是不稳定的算法对问题1进行计算,并与真解作比较。(选改进的欧拉法)2.1实验思路:算法的稳定性是与步长h密切相关的。而对于问题一而言,取定步长h=0.1不论是单步法或低阶多步法都是稳定的算法。所以考虑改变h取值范围,借此分析不同步长会对结果造成什么影响。故依次采用h=2.0、2.2、2.4、2.6的改进欧拉法。2.2实验代码:%主程序x=0:3:30;y=exp(-x)+x;dyfun=inline(-y+x+1); x1,m1,y1=naeuler2(dyfun,0,20,1,2);x2,m2,y2=naeuler2(dyfun,0,22,1,2.2);x3,m3,y3=naeuler2(dyfun,0,24,1,2.4);x4,m4,y4=naeuler2(dyfun,0,26,1,2.6);subplot(2,2,1)plot(x,y,r,x1,y1,+);xlabel(h=2.0);subplot(2,2,2)plot(x,y,r,x2,y2,+);xlabel(h=2.2);subplot(2,2,3)plot(x,y,r,x3,y3,+);xlabel(h=2.4);subplot(2,2,4)plot(x,y,r,x4,y4,+);xlabel(h=2.6);2.3实验结果:xh=2.0h=2.2h=2.4h=2.60.0 1.0000 1.0000 1.0000 1.0000 0.1 3.0000 3.4200 3.8800 4.3800 0.2 5.0000 5.8884 6.9904 8.3684 0.3 7.0000 8.4158 10.4418 13.4398 0.4 9.0000 11.0153 14.3979 20.4388 0.5 11.0000 13.7027 19.1008 30.8690 0.6 13.0000 16.4973 24.9092 47.4068 0.7 15.0000 19.4227 32.3536 74.8161 0.8 17.0000 22.5077 42.2194 121.5767 0.9 19.0000 25.7874 55.6687 202.7825 1.0 21.0000 29.3046 74.4217 345.3008 实验结果分析:从实验1结果可以看出,在算法满足收敛性和稳定性的前提下,Eluer法虽然计算并不复杂,凡是精度不足,反观改进的Eluer法和SK法虽然计算略微复杂但是结果很精确。实验2改变了步长,导致算法理论上收敛但是不满足稳定性。结果表示步长h越大,结果越失真。对于同一个问题,步长h的选取变得尤为重要,这三种单步法的绝对稳定区间并不一样,所以并没有一种方法是万能的,我们应该根据不同的步长来选取合适的方法。实验二:Ritz-Galerkin方法与有限差分法1、 用中心差分格式求解边值问题取步长h=0.1,并与真解作比较。1.1实验代码:%中心差分法function U=fdm(xspan,y0,y1,h)%xspan为x取值范围,y0,y1为边界条件,h为步长N=1/h;d=zeros(1,N-1);for i=1:N x(i)=xspan(1)+i*h; q(i)=1; f(i)=x(i);endfor i=1:N-1 d(i)=q(i)*h*h+2;end a=diag(d); b=zeros(N-1); c=zeros(N-1);for i=1:N-2 b(i+1,i)=-1;endfor i=1:N-2 c(i,i+1)=-1;endA=a+b+c;for i=2:N-2 B(i,1)=f(i)*h*h;end B(1,1)=f(1)*h*h+y0; B(N-1,1)=f(N-1)*h*h+y1; U= inv(A)*B;%主程序x=0:0.1:1;y=x+(exp(1)*exp(-x)/(exp(2)-1)-(exp(1)*exp(x)/(exp(2)-1);y1=fdm(0,1,0,0,0.1);y1=0,y1,0;plot(x,y,r,x,y1,+)xlabel(x);ylabel(y);legend(y真解,y1中心差分法,Location,NorthWest);1.2实验结果:xy真解y1中心差分法0.0 0.0000 0.0000 0.1 0.0148 0.0148 0.2 0.0287 0.0287 0.3 0.0409 0.0408 0.4 0.0505 0.0504 0.5 0.0566 0.0565 0.6 0.0583 0.0582 0.7 0.0545 0.0545 0.8 0.0443 0.0443 0.9 0.0265 0.0265 1.0 0.0000 0.0000 2、用Ritz-Galerkin方法求解上述问题,并与真值作比较,列表画图。2.1实验代码:%Ritz_Galerkin法function vu=Ritz_Galerkin(x0,y0,x1,y1,h)%x0,x1为x取值范围,y0,y1为边界条件,h为步长N=1/h;syms x;for i=1:N fai(i)=x*(1-x)*(x(i-1); dfai(i)=diff(x*(1-x)*(x(i-1); endfor i=1:N for j=1:N fun=dfai(i)*dfai(j)+fai(i)*fai(j); A(i,j)=int(fun,x,0,1); end fun=x*fai(i)+dfai(i); f(i)=int(fun,x,0,1);endc=inv(A)*f;product=c.*fai; sum=0; for i=1:N sum=sum+product(i);endvu=;for y=0:h:1 v=subs(sum,x,y); vu=vu,v; endy=0:h:1;yy=0:0.1:1; u=sin(yy)/sin(1)-yy; u=vpa(u,5);vu=vpa(vu,5); %主程序x=0:0.1:1;y=x+(exp(1)*exp(-x)/(exp(2)-1)-(exp(1)*exp(x)/(exp(2)-1);y1=Ritz_Galerkin(0,0,1,0,0.1);y1=double(y1);plot(x,y,r,x,y1,+)xlabel(x);ylabel(y);legend(y为真解,y1为RG法,Location,NorthWest);2.2实验结果:xy真解y1RG法0.0 0.0000 0.0000 0.1 0.0148 0.0148 0.2 0.0287 0.0287 0.3 0.0409 0.0409 0.4 0.0505 0.0505 0.5 0.0566 0.0566 0.6 0.0583 0.0583 0.7 0.0545 0.0545 0.8 0.0443 0.0443 0.9 0.0265 0.0265 1.0 0.0000 0.0000 3、若
温馨提示
- 1. 本站所有资源如无特殊说明,都需要本地电脑安装OFFICE2007和PDF阅读器。图纸软件为CAD,CAXA,PROE,UG,SolidWorks等.压缩文件请下载最新的WinRAR软件解压。
- 2. 本站的文档不包含任何第三方提供的附件图纸等,如果需要附件,请联系上传者。文件的所有权益归上传用户所有。
- 3. 本站RAR压缩包中若带图纸,网页内容里面会有图纸预览,若没有图纸预览就没有图纸。
- 4. 未经权益所有人同意不得将文件中的内容挪作商业或盈利用途。
- 5. 人人文库网仅提供信息存储空间,仅对用户上传内容的表现方式做保护处理,对用户上传分享的文档内容本身不做任何修改或编辑,并不能对任何下载内容负责。
- 6. 下载文件中如有侵权或不适当内容,请与我们联系,我们立即纠正。
- 7. 本站不保证下载资源的准确性、安全性和完整性, 同时也不承担用户因使用这些下载资源对自己和他人造成任何形式的伤害或损失。
最新文档
- 医生个人先进事迹
- 宿州市中石油2025秋招面试半结构化模拟题及答案机械与动力工程岗
- 厨房承包合同集锦15篇
- 2025年潍坊诸城市市属国有企业公开招聘工作人员(9名)模拟试卷完整参考答案详解
- 2025年甘肃省天水天光半导体有限责任公司招聘18人模拟试卷完整参考答案详解
- 2025年山东省药品不良反应监测中心公开招聘人员考前自测高频考点模拟试题及答案详解(历年真题)
- 2025年甘肃省大数据中心招聘工作人员模拟试卷及答案详解(网校专用)
- 2025年临沂市罗庄区教育系统部分事业单位公开招聘教师(43名)模拟试卷附答案详解(模拟题)
- 2025年上半年四川乐山职业技术学院赴四川大学考核招聘10人模拟试卷完整答案详解
- 2025年武汉工程大学人才引进33人模拟试卷及答案详解(各地真题)
- 陕西省西工大附中2022-2023学年七年级上学期第一次月考英语试卷(含答案)
- 2025辅警考试题库(含答案)
- QGDW10212-2019电力系统无功补偿技术导则
- 牛奶面包食品配送服务 投标方案(技术方案)
- 菜鸟驿站运营管理制度
- 风管安装协议书
- 个人车位租赁合同(含充电桩安装)
- 2025年人教版小学六年级上册奥林匹克数学竞赛测试题(附参考答案)
- 第三方垫付协议合同协议
- 《鸿蒙HarmonyOS应用开发基础》全套教学课件
- 订购包装木箱合同协议
评论
0/150
提交评论