(精品word)MATLAB龙格-库塔方法解微分方程.doc

(精品word)MATLAB龙格-库塔方法解微分方程
(精品word)MATLAB龙格-库塔方法解微分方程
(精品word)MATLAB龙格-库塔方法解微分方程
龙格—库塔方法是一种经典方法，具有很高的精度，它间接的利用了泰勒级数展开，避免了高阶偏导数的计算。此处以最为经典的四级四阶龙格-库塔方法为例,计算格式如下

1龙格—库塔法解一阶ODE
对于形如的一阶ODE初值问题，可以直接套用公式，如今可以借助计算机方便的进行计算，下面给出一个实例

取步长h=0。1,此处由数学知识可得该方程的精确解为。在这里利用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)编写主函数解微分方程,并观察数值解与精确解的差异

clearall
h=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;
y(k+1）=runge_kutta（h，x（k），y(k)）；
end
z=sqrt(1+2＊x)；
plot（x，y,’*’)；
holdon
plot（x，z,'r'）;

结果如下图,数值解与解析解高度一致

2龙格—库塔法解高阶ODE

对于高阶ODE来说，通用的方法是将高阶方程通过引入新的变量降阶为一阶方程组，此处仍以一个实例进行说明。

这是一个二阶ODE，描述的是一个物体的有阻尼振动情况。
初始条件为,将方程降阶，引入一个向量型变量Y
故有
记则至此，二阶方程降阶为一阶方程组。值得注意的是此时再用龙格-库塔法进行求解时，代入的将是一个Y向量。同样利用MATLAB进行计算,步长h=0。05,时间周期为［0,20].

编写ODE函数
functionY=odefun1（~,Y0)
%此处Y0为一个列向量，因为时间t未显含在一阶方程组中
％所以ode函数的第一个参数为空，要根据具体情况而定。
Y=［Y0(2）；
(2000-200＊Y0（2）-750*Y0（1)）/500;]；
end

编写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＊k2）;
k4=odefun1(t0+h,Y0+h＊k3)；
Y=Y0+h＊（k1+2＊k2+2＊k3+k4)/6;
end

编写主函数

clearall
h=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）；

end
plot（t,z(1,：)，'r'）；％位移y的图像
holdon
plot（t，z（2,：)）；%速度y’的图像

求解结果如下图