寸_铁

这个屌丝很懒，什么也没留下！

热门标签

matlab求解常微分方程初值问题的三类方法_常微分方程初值问题matlab

作者：寸_铁 | 2024-07-12 12:45:08

踩

常微分方程初值问题matlab

求微分方程初值问题： $y''-(1-y^2)y'+y=0,y(0)=2,y'(0)=0$

我们以二阶微分方程为例来介绍这三种方法的使用，当然，第一种方法——四阶龙格-库塔法实际上就是自己手动编写程序来求解，而第二种方法——ode45就是matlab上的也是用四阶龙格-库塔法实现的更完备的方法，都是数值解法，第三种方法是微分方程的符号解法，是精确解法。

讨论一阶微分方程

$\frac{dy}{dt}=f(t,y)$

初始条件为 $y(t_0)=y_0$ ，其中y和f可以是矢量（即微分方程组）。

所谓的数值积分，是解决在初值已知的情况下，对 $f(t,y)$ 进行近似积分的问题，通常称为微分方程的初值问题。

对上述一阶微分方程两边进行积分：

$y(t)=y(t_0)+\int_{t_0}^t {f(\tau,y)} \,{\rm d}\tau$

在离散的时间下： $t=t_0,t_1,t_2,\cdots,t_m,t_{m+1}$

$y(t_{m+1})=y(t_0)+\int_{t_0}^{t_{m+1}} {f(\tau,y)} \,{\rm d}\tau=y(t_m)+\int_{t_m}^{t_{m+1}} {f(\tau,y)} \,{\rm d}\tau$

注：如果要求解的是高阶微分方程，可降阶为一阶微分方程组。如二阶微分方程 $y''(t)=f(t,y)$

可降阶为：

$y_1'=y_2$

$y_2'=f(t,y_1)$

也就是把原来二阶的 $y$ 看成 $y_1$ ： $y_1''(t)=f(t,y_1)$ ，将 $y_1$ 一阶导数令为 $y_2$ ： $y_2'=f(t,y_1)$

步长： $h=t_{m+1}-t_m$

单步法：只由前一时刻的数值 $y(t_m)$ 就可以求得后一时刻的数值 $y(t_{m+1})$ ，是一种自启动算法。（如欧拉法、龙格-库塔法等）

多步法：计算 $y(t_{m+1})$ 时要用到 $t_m,t_{m-1},t_{m-2},\cdots$ y的数值，这种方法为多步法，不是自启动算法。（如Adams法）

截断误差：分析数值积分的精度常用泰勒级数，设 $y_m$ 准确，则 $y_{m+1}$ :

$y(t_{m+1})=y(t_m)+y'(t_m)h+1/2! y''(t_m)h^2+\cdots+1/r!y^{(r)}(t_m)h^r+o(h^{r+1})$

若只取前几项之和来近似计算 $y_{m+1}$ ，这样产生的误差称为截断误差，比如取到 $1/k! y^{(k)}(t_m)h^k$ ，就称具有k阶精度，即该方法是k阶的。

一、四阶龙格-库塔法

在看四阶之前，先了解一下二阶的龙格-库塔法。

二阶，即泰勒展开到二阶：

$y(t_{m+1})=y(t_m)+hy'(t_m)+\frac{h^2}{2}y''(t_m)$

将

$y'(t_m)=f(t_m,y_m),y''(t_m)=(\frac{\partial f}{\partial t} +f\frac{\partial f}{\partial y})|_{t=t_m}$

代入上式，然后假设微分方程的解可写成如下形式：

$y_{m+1}=y_m+(a_1k_1+a_2k_2)h$

$k_1=f(t_m,y_m)$

$k_2=f(t_m+b_1h,y_m+b_2k_1h)$

取不同的 $a_1,a_2,b_1,b_2$ 可以得到不同的二阶龙格-库塔法。

那么四阶龙格-库塔法类似：

$y_{m+1}=y_m+\frac{h}{6}(k_1+2k_2+2k_3+k_4)$

$k_1=f(t_m,y_m)$

$k_2=f(t_m+h/2,y_m+k_1h/2)$

$k_3=f(t_m+h/2,y_m+k_2h/2)$

$k_4=f(t_m+h,y_m+k_3h)$

具体计算时，由 $y_0$ 可以求得 $y_1$ ，由 $y_1$ 可以求得 $y_2$ ，……，直到所要求的值。

附上matlab代码框架：

首先定义一阶微分方程组：


function func=f(t,y)
f1=y(2);
f2=(1-y(1)^2)*y(2)-y(1);
func=[f1;f2];
end

然后是龙格-库塔算法：


function [t,Y]=RKutta(Delta,Y0,tm)
% Delta为步长，Y0为初始值，是一个列向量,tm是计算
Y(:,1)=Y0;
k=1;
t=0:Delta:tm;
t_size=length(t);
for i=1:t_size
    z1=f(t(k),Y(:,k));
    z2=f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
    z3=f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
    z4=f(t(k)+Delta,Y(:,k)+z3*Delta);
    Y(:,k+1)=Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
    k=k+1;
end
Y=Y(:,1:end-1); % 其实求得的是0到20+Delta秒的y
end

设定初值和步长，进行计算：


Delta=0.001;
Y0=[2;0];
tm=20;
[t1,y1]=RKutta(Delta,Y0,tm)

将y随时间变化规律画出来：

二、ode45

matlab提供了好几种求解器，其中最常用的就是ode45求解器。

其最简单的使用方法是：[t,y] = ode45(odefun,tspan,y0)

同样首先定义求解的微分方程组odefun：


function dydt = vdp1(t,y)
 
dydt = [y(2); (1-y(1)^2)*y(2)-y(1)];
 
end

求解时间范围tspan=[0,20]，初值y0=[2;0]

[t,y] = ode45(@vdp1,[0 20],[2; 0]);

求得的是y随着t变化的序列。

同样画出图像：

此外，如果积分区间不定，但知道需要在y减小到-1.5时停止积分，就需要使用ode45的事件描述：


function [position,isterminal,direction] = myEventsFcn(t,y)
position = y(1)+1.5; % 该式等于0即事件发生（即y(1)=-1.5）
isterminal = 1;  % 当事件多于1个时，设置当某个事件发生时停止积分（为1）
direction = -1;   % 值为 -1 时仅在事件函数递减的位置查找零

ode45的事件设置：


options=odeset('events',@myEventsFcn);
[t2,y2] = ode45(@vdp1,[0 20],[2; 0],options);

此时画出图像为：

三、dsolve

符号解也是精确解，然而大多数方程并不能求得解析解，比如本例中的这个二阶微分方程，用此方法求解会得到报错：


syms t y(t)
eqn=diff(y,t,2)-(1-y^2)*diff(y,t,1)+y==0;
Dy = diff(y,t);
cond=[y(0)==2,Dy(0)==0];
ySol(t)= dsolve(eqn,cond)

所以还是需要用数值方法去求解。

这里举一个可以求解精确解的例子：


syms t y(t)
eqn=diff(y,t,2)-3*diff(y,t,1)+2*y==1;
Dy = diff(y,t);
cond=[y(0)==1,Dy(0)==0];
ySol(t)= dsolve(eqn,cond)

声明：本文内容由网友自发贡献，不代表【wpsshop博客】立场，版权归原作者所有，本站不承担相应法律责任。如您发现有侵权的内容，请联系我们。转载请注明出处：https://www.wpsshop.cn/w/寸_铁/article/detail/814155