MATLAB龙格-库塔方法解微分方程.docx

龙格-库塔方法是一种经典方法,具有很高的精度,它间接的利用了泰勒级数展开,避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例,计算格式如下1龙格-库塔法解一阶ODE对于形如的一阶ODE初值问题,可以直接套用公式,如今可以借助计算机方便的进行计算,下面给出一个实例取步长h=,此处由数学知识可得该方程的精确解为。在这里利用MATLAB编程,计算数值解并与精确解相比,代码如下:(1)写出微分方程,便于调用和修改functionval=odefun(x,y)val=y-2*x/y;end(2)编写runge-kutta方法的函数代码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=;x0=0;y0=1;x=:h:1;y(1)=runge_kutta(h,x0,y0);fork=1:length(x)x(k)=x0+k*h;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向量。同样利用MATLAB进行计算,步长h=,时间周期为[0,20].编写ODE函数functionY=odefun1(~,Y0)%此处Y0为一个列向量,因为时间t未显含在一阶方程组中%所以ode函数的第一个参数为空,要根据具体情况而定。Y=[Y0(2);(2000-200*Y0(2)-750*Y0(1))/500;