资源描述：
《MATLAB龙格-库塔方法解微分方程.docx》由会员上传分享，免费在线阅读，更多相关内容在工程资料-天天文库。
1、龙格-库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例，计算格式如下1龙格-库塔法解一阶ODE对于形如的一阶ODE初值问题，可以直接套用公式，如今可以借助计算机方便的进行计算，下面给出一个实例取步长h=0.1，此处由数学知识可得该方程的精确解为。在这里利用MATLAB编程，计算数值解并与精确解相比，代码如下：(1)写出微分方程，便于调用和修改functionval=odefun(x,y)val=y-2*x/y;end(2)编写runge-kutta方法的函数
2、代码functiony=runge_kutta(h,x0,y0)k1=odefun(x0,y0);k2=odefun(x0+h/2,y0+h/2*k1);k3=odefun(x0+h/2,y0+h/2*k2);k4=odefun(x0+h,y0+h*k3);y=y0+h*(k1+2*k2+2*k3+k4)/6;end(3)编写主函数解微分方程,并观察数值解与精确解的差异clearallh=0.1;x0=0;y0=1;x=0.1:h:1;y(1)=runge_kutta(h,x0,y0);fork=1:length(x)x(k)=x0+k*h;
3、y(k+1)=runge_kutta(h,x(k),y(k));endz=sqrt(1+2*x);plot(x,y,’*’);holdonplot(x,z,'r');结果如下图，数值解与解析解高度一致2龙格-库塔法解高阶ODE对于高阶ODE来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处仍以一个实例进行说明。这是一个二阶ODE,描述的是一个物体的有阻尼振动情况。初始条件为，将方程降阶，引入一个向量型变量Y故有记则至此，二阶方程降阶为一阶方程组。值得注意的是此时再用龙格-库塔法进行求解时，代入的将是一个Y向量。同样利用MATL
4、AB进行计算，步长h=0.05，时间周期为[0,20].(1)编写ODE函数functionY=odefun1(~,Y0)%此处Y0为一个列向量，因为时间t未显含在一阶方程组中%所以ode函数的第一个参数为空，要根据具体情况而定。Y=[Y0(2);(2000-200*Y0(2)-750*Y0(1))/500;];end(2)编写runge-kutta函数functionY=rkfa(h,t0,Y0)k1=odefun1(t0,Y0);k2=odefun1(t0+h/2,Y0+h/2*k1);k3=odefun1(t0+h/2,Y0+h/2*k
5、2);k4=odefun1(t0+h,Y0+h*k3);Y=Y0+h*(k1+2*k2+2*k3+k4)/6;end(1)编写主函数clearallh=0.05;t=0.05:h:20;t0=0;Y0=[0;0];%初值Y=cell(1,length(t));Y{1}=rkfa(h,t0,Y0);z=zeros(2,length(t));fork=1:length(t)Y{k+1}=rkfa(h,t0,Y{k});z(1,k)=Y{k}(1);z(2,k)=Y{k}(2);endplot(t,z(1,:),'r');%位移y的图像holdon
6、plot(t,z(2,:));%速度y’的图像求解结果如下图
