matlab 四阶龙格-库塔法求微分方程

Matlab 实现四阶龙格-库塔发求解微分方程

从理论上讲，只要函数在某区间上充分光滑，那么它可以展开为泰勒级数，因此在该区间上的函数值可用各阶导数值近似地表示出来，反之其各阶导数值也可用某些函数值的线性组合近似地表示出来。龙格-库塔法就是将待求函数y(t)展开为泰勒级数，并用方程函数f(t,y)近似其各阶导数，从而迭代得到y(t)的数值解。具体来说，四阶龙格-库塔迭代公式为

1yn?1?yn?h(k1?2k2?2k3?k4)

6k1?f(tn,yn)

k2?f(tn?h/2,yn?hk1/2) k3?f(tn?h/2,yn?hk2/2) k3?f(tn?h,yn?hk3)

实验内容：

?1?x2，x?2??0.4x1?0.2x2?0.5u，x1(0)?x2(0)?0，u为单位阶已知二阶系统x跃信号。用四阶龙格-库塔法求数值解。分析步长对结果的影响。 实验总结：

实验报告要求简要的说明实验原理；简明扼要地总结实验内容；编制m文件，并给出运行结果。报告格式请按实验报告模板编写。

进入matlab，

Step1：choose way1 or way2

way1）：

可以选择直接加载M文件（函数M文件）。

way2）：

点击new——function，先将shier（函数1文本文件）复制运行; 点击new——function，再将RK（函数2文本文件）运行; 点击new——function，再将finiRK（函数3文本文件）运行;

Step2：

回到command页面输入下面四句。

[t,k]=finiRK45([0;0],150);%迭代150次，步长=20/150

[t1 k1]=ode45(@shier,[0 -10],[0 0]);%调用matlab自带四阶龙格-库塔，对比结果

[t2 k2]=ode45(@shier,[0 10],[0 0]);

plot(t,k(1,:),'-',t1,k1(:,1),'*',t2,k2(:,1),'^')%在图形上表示出来

补充：改变步长影响数据的准确性。

函数1 shier：

function dx =shier(t,x)

%UNTITLED Summary of this function goes here % Detailed explanation goes here dx=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为上限（为了方便fniRK45这 %里默认下限为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)üh为迭代初值，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 1] y11=fliplr(y11);

x1=[t11 t12];%将所有的解和成最终的解向量 y1=[y11 y12]; end

%步长=(10+10)/N1

本文来源：https://www.bwwdw.com/article/v5pf.html