Matlab求常微分方程组的解析解

MATLAB代码 专栏收录该内容
21 篇文章 4 订阅

最近同学毕设需要求解循坏摆的微分方程,我在帮忙过程中学习了一下常微分方程的解析解和数值解的求法,在此分享。
以下讲解遵循Matlab官方文档提供的方程和写法。
(强烈建议大家有问题多看官方文档,非常有用)

介绍一下核心函数

常用形式:S = dsolve(eqn,cond)
功能介绍:解微分方程eqn,其中eqn是一个符号方程,cond是初值条件(边界条件)。使用diff和==表示微分方程。例如,diff(y,x)==y表示方程dy/dx=y。通过将eqn指定为这些方程的向量来求解微分方程组。

1. 一阶常微分方程求解(无初值)

方程:dy/dt=ay
代码

syms y(t) a; %syms 定义方程用的变量和未知常量(可有可无)。y(t)表示y是关于t的函数,求解t。 
eqn=diff(y,t)==a*y; %diff(y,t)即dy/dt,默认为1阶。
S=dsolve(eqn) %dsolve()为求解函数,将方程放入即可

说明:首先定义方程所用变量和未知常量,然后列出微分方程,最后dsolve()求解。
结果

S =
C2*exp(a*t)

分析:C2为未定的常数项,需要提供一个初值,才能确定。

2. 二阶常微分方程求解(无初值)

方程:dy²/dt²=ay
代码

syms y(t) a; %syms 定义方程用的变量和未知常量。y(t)表示y是关于t的函数,求解t。 
eqn=diff(y,t,2)==a*y; %diff(y,t,2)即dy²/dt²,第三个参数表示阶数,更高阶数一样。
S=dsolve(eqn) %dsolve()为求解函数,将方程放入即可

结果

S =
C3*exp(a^(1/2)*t) + C4*exp(-a^(1/2)*t)

分析:无初值情况下,C3,C4未定,更高阶数解法类推,修改diff()第三个参数即可。

3. 二阶微分方程求解(有初值)

方程:dy²/dt²=ay
代码

syms y(t) a; %syms 定义方程用的变量和未知常量。y(t)表示y是关于t的函数,求解t。 
eqn=diff(y,t,2)==a*y; %定义diff(y,t,2)即dy²/dt²,第三个参数表示阶数。
Dy=diff(y,t); %定义方程Dy为dy/dt;
cond=[y(0)==5,Dy(0)==1]; %定义初值,t=0时,y=5,dy/dt=1;
S=dsolve(eqn,cond) %dsolve()为求解函数,将方程放入即可

说明:一阶微分方程求解时,不需要定义Dy,二阶方程有dy/dt的初值时,需要定义
结果

S =
(exp(a^(1/2)*t)*(5*a^(1/2) + 1))/(2*a^(1/2)) + (exp(-a^(1/2)*t)*(5*a^(1/2) - 1))/(2*a^(1/2))

分析:有初值情况下,无未定常数项

4. 一阶微分方程组求解(有初值)

方程:dy/dt=z;
dz/dt=-y
代码

syms y(t) z(t); %syms 定义方程用的变量 
eqn=[diff(y,t)==z,diff(z,t)==-y]; %定义方程组dy/dt=z,dz/dt=-y
cond=[y(0)==1,z(1)==1]; %定义初值,t=0时,y=5;t=1时,y=1;
[ySol zSol]=dsolve(eqn,cond) %dsolve()为求解函数,将方程放入即可

解释:老套路,只不过要定义俩方程,给俩初值
结果

ySol =
cos(t - 1)/(2*cos(1)) - cos(t + 1)/(2*cos(1)) + cos(t - atan(1/cos(1)))*(1/cos(1)^2 + 1)^(1/2)
zSol =
sin(t + 1)/(2*cos(1)) - sin(t - 1)/(2*cos(1)) + (cos(t + atan(cos(1)))*(cos(1)^2 + 1)^(1/2))/cos(1)

5. 无法求解情况

方程:dy²/dt²=(1-y²)*dy/dt-y
代码

syms y(t) %syms 定义方程用的变量和常量
eqn=diff(y,2)==(1-y^2)*diff(y)-y; %定义方程dy²/dt²=(1-y²)*dy/dt-y
cond=y(0)==1; %定义初值,t=0时,y=1
S=dsolve(eqn,cond) %dsolve()为求解函数,将方程放入即可

解释:方程本身看起来不复杂,但是无法求解
结果

警告: Explicit solution could not be found. 
> In dsolve (line 201) 
S =
[ empty sym ]

分析:特定解无法求出,结果是空的方程

6. 无法求解解决方案

遇到无法求解的情况,dsolve就束手无策了,但是Matlab依旧强大,提供了求数值解的方法,非常有效。
我会在另一篇博文里面讲解如何对复杂的微分方程进行数值求解,这才是我解决同学毕设问题的方案。(已写,点击此处查看)
这里先把同学的题目放在这里,大家可以编写解析解求解方案试一下,反正我求解不出来。
循环摆微分方程如下:
在这里插入图片描述
给定的初值(w接近0,但实际上不能设置为0):
在这里插入图片描述
我的求解代码如下:

syms R(t) v(t) w(t);
eqns=[diff(R,t)==-v/w-0.01,...
    diff(v,t)==(10*9.8-(w^2*R+9.8*sin(t))*exp(0.1*(t+pi/2)))/(w*11),...
    diff(w,t)==w*0.01/R+2*v/R+9.8*cos(t)/(w*R)];
cond=[R(0)==1,v(0)==0,w(0)==1e3];
S=dsolve(eqns,cond);

7. 加油!

  • 9
    点赞
  • 0
    评论
  • 59
    收藏
  • 打赏
    打赏
  • 扫一扫,分享海报

参与评论 您还未登录,请先 登录 后发表或查看评论
©️2022 CSDN 皮肤主题:大白 设计师:CSDN官方博客 返回首页

打赏作者

y=520(2sinM+sin2M)

你的鼓励将是我创作的最大动力

¥2 ¥4 ¥6 ¥10 ¥20
输入1-500的整数
余额支付 (余额:-- )
扫码支付
扫码支付:¥2
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值