四阶龙格库塔法求解常微分方程的初值问题matlab通用程序.docx

下载文档 降价啦

15
0
约2.75千字
约 3页
2017-07-19 发布于浙江
举报
版权申诉
保障服务

四阶龙格库塔法求解常微分方程的初值问题matlab通用程序.docx

1、本文档共3页，可阅读全部内容。
2、有哪些信誉好的足球投注网站（book118）网站文档一经付费（服务费），不意味着购买了该文档的版权，仅供个人/单位学习、研究之用，不得用于商业用途，未经授权，严禁复制、发行、汇编、翻译或者网络传播等，侵权必究。
3、本站所有内容均由合作方或网友上传，本站不对文档的完整性、权威性及其观点立场正确性做任何保证或承诺！文档内容仅供研究参考，付费前请自行鉴别。如您付费，意味着您自己接受本站规则且自行承担风险，本站不退款、不进行额外附加服务；查看《如何避免下载的几个坑》。如果您已付费下载过本站文档，您可以点击这里二次下载。
4、如文档侵犯商业秘密、侵犯著作权、侵犯人身权等，请点击“版权申诉”（推荐），也可以打举报电话：400-050-0827(电话支持时间：9:00-18:30)。

四阶龙格库塔法求解常微分方程的初值问题matlab通用程序

参考教材《数值分析》李乃成.梅立泉clearclcformat longm=input(请输入常微分方程的阶数m=);a=input(请输入x下限a=);b=input(请输入x上限b=);h=input(请输入步长h=);ym=input(令y(1,1)=y,y（2,1）=y’,y(3,1)=y’’...请输入ym=,s); ? ? %输入的时候必须按照这个形式输入y1=y(1,1);if m==1 ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %一阶初值问题单独求解? ? mm=(b-a)/h;? ? y(1,1)=input(请输入在初值点的函数值f(a)=);? ? x=a;? ? y11(1)=y(1,1);? ? for k1=2:(mm+1)? ? ? ? y1=y(1,1);? ? ? ? K(1,1)=h*(eval(ym)); ? ? ? ? ? ? ? ? ? ? ? ? %计算K1? ? ? ? x=x+h/2;? ? ? ? y(1,1)=y1+K(1,1)/2;? ? ? ? y1=y(1,1);? ? ? ? K(1,2)=h*(eval(ym)); ? ? ? ? ? ? ? ? ? ? ? ? %计算K2? ? ? ? x=x;? ? ? ? y(1,1)=y1+K(1,2)/2-K(1,1)/2;? ? ? ? y1=y(1,1);? ? ? ? K(1,3)=h*(eval(ym)); ? ? ? ? ? ? ? ? ? ? ? ? ?%计算K3? ? ? ? x=x+h/2;? ? ? ? y(1,1)=y1+K(1,3)-K(1,2)/2;? ? ? ? y1=y(1,1);? ? ? ? K(1,4)=h*(eval(ym)); ? ? ? ? ? ? ? ? ? ? ? ? ?%计算K4? ? ? ? y11(k1)=y11(k1-1)+(K(1,1)+2*K(1,2)+2*K(1,3)+K(1,4))/6;? ? ? ? y(1,1)=y11(k1);? ? ? ? x=a+(k1-1)*h;? ? ? ?? ? endy11else ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %高阶初值问题? ? mm=(b-a)/h; ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %一共要求解mm个数据点? ? for k2=1:m ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?%读取初值条件? ? ? ? fprintf(请输入%d阶导数的初值f(%d)(a)=\n,(k2-1),(k2-1));? ? ? ? y(k2,1)=input(=);? ? end? ? for k2=1:m ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ?y22(1,k2)=y(k2,1); ? ? ? ? ? ? ? ? ? ? ? ? ?%先把初值保存在矩阵y22(m,n)中,m表示第几个所求点，n表示第n阶初值? ? end? ? x=a;? ? for k4=2:(mm+1) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %求解mm个数据点的循环? ? ? ? for k=1:(m-1) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? %计算K1，包括每一阶的K1? ? ? ? ? ? K(k,1)=h*y(k+1,1); ? ? ? ? ? ? ? ? ? ? ? ?%y(k+1,1)中k+1表示第k+1阶，1表示第一个点；K(k,1)中k表示阶数，1表示K1? ? ? ? end? ? ? ? K(m,1)=h*(eval(ym));? ? ? ? x=x+h/2; ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?%求解K1之前，先重新对x和y赋值? ? ? ? for k3=1:m ? ? ? ? ? ? ?? ? ? ? ? ? y(k3,1)=y(k3,1)+K(k3,1)/2;? ? ? ? end? ? ? ? for k=1:(m-1) ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ?%计算K2? ? ? ? ? ? K(k,2)=h*y(k+1,1);? ? ? ? end? ? ? ? K(m,2)=h*(eval(ym));? ? ? ? x=x;? ? ? ? for k3=1:m? ? ? ? ? ? y(k3,1)=y(k3,1)-K(k3,1)/2+K(k3,2)/2;? ? ? ? end?