《matlab 四阶龙格-库塔法求微分方程(5页).doc》由会员分享,可在线阅读,更多相关《matlab 四阶龙格-库塔法求微分方程(5页).doc(5页珍藏版)》请在taowenge.com淘文阁网|工程机械CAD图纸|机械工程制图|CAD装配图下载|SolidWorks_CaTia_CAD_UG_PROE_设计图分享下载上搜索。
1、- Matlab 实现四阶龙格-库塔发求解微分方程从理论上讲,只要函数在某区间上充分光滑,那么它可以展开为泰勒级数,因此在该区间上的函数值可用各阶导数值近似地表示出来,反之其各阶导数值也可用某些函数值的线性组合近似地表示出来。龙格-库塔法就是将待求函数展开为泰勒级数,并用方程函数近似其各阶导数,从而迭代得到的数值解。具体来说,四阶龙格-库塔迭代公式为实验内容:已知二阶系统,u为单位阶跃信号。用四阶龙格-库塔法求数值解。分析步长对结果的影响。实验总结:实验报告要求简要的说明实验原理;简明扼要地总结实验内容;编制m文件,并给出运行结果。报告格式请按实验报告模板编写。进入matlab,Step1:c
2、hoose way1 or way2way1):可以选择直接加载M文件(函数M文件)。way2):点击newfunction,先将shier(函数1文本文件)复制运行;点击newfunction,再将RK(函数2文本文件)运行;点击newfunction,再将finiRK(函数3文本文件)运行;Step2:回到command页面输入下面四句。t,k=finiRK45(0;0,150);%迭代150次,步长=20/150t1 k1=ode45(shier,0 -10,0 0);%调用matlab自带四阶龙格-库塔,对比结果t2 k2=ode45(shier,0 10,0 0);plot(t,k(
3、1,:),-,t1,k1(:,1),*,t2,k2(:,1),)%在图形上表示出来补充:改变步长影响数据的准确性。函数1 shier:function dx =shier(t,x)%UNTITLED Summary of this function goes here% Detailed explanation goes heredx=zeros(2,1);dx(1)=x(2);dx(2)=-0.4*x(1)-0.2*x(2)+0.5*0.5*(sign(t)-sign(-t);end函数2 RK45:function t,y=RK45(f,b,ch0,N)%f为函数句柄,b为上限(为了方便f
4、niRK45这%里默认下限为0),ch0为初值,N为迭代次数。h=b/N;%计算步长y=zeros(2,N);%开辟y,k的向量空间确定维数t=zeros(1,N);t(1)=0;y(:,1)=ch0;%赋迭代初值for j=1:N %循环迭代过程 t(j+1)=t(j)+h ; k1=f(t(j),y(:,j); k2=f(t(j)+h/2,y(:,j)+h*k1/2); k3=f(t(j)+h/2,y(:,j)+h*k2/2);k4=f(t(j)+h/2,y(:,j)+h*k3/2); y(:,j+1)=y(:,j)+h*(k1+2*k2+2*k3+k4)/6;end end函数3 finiRK45:function x1, y1 =finiRK45(fch0,N1)%fch为迭代初值,N1为迭代次数t11,y11=RK45(shier,-10,fch0,N1/2);%求在【-10,10】的解t12,y12=RK45(shier,10,fch0,N1/2);t11=fliplr(t11);%左右转换如:a1=1 2 3执行后a1=3 2 1y11=fliplr(y11);x1=t11 t12;%将所有的解和成最终的解向量y1=y11 y12;end%步长=(10+10)/N1第 5 页-