• 设为首页
  • 点击收藏
  • 手机版
    手机扫一扫访问
    迪恩网络手机版
  • 关注官方公众号
    微信扫一扫关注
    公众号

用Matlab求解微分方程

原作者: [db:作者] 来自: [db:来源] 收藏 邀请

用Matlab求解微分方程

解微分方程有两种解,一种是解析解,一种是数值解,这两种分别对应不同的解法

解析解

利用dsolve函数进行求解

syms  x;
s = dsolve(\'eq1,eq2,...\', ’cond1,cond2,...\', \'v\');
%eq:微分方程
%cond:条件
%v:独立变量
%形如:方程:y\'= f(t,y),初值:y(t0) = y0

1.求解析解

 dsolve(\'Du = 1+ u^2\',\'t\')
 ans =
 
 tan(C2 + t)
          1i
         -1i

求 的解析解

s = dsolve(\'D2y=3*y+2*x\',\'x\');
% D2y用以表示y的二阶导数,默认是以t为自变量的,所以最好指明自变量为x.
syms y(x);
s = dsolve([diff(y,x,2) == 3*y+2*x], [y(0) == 5]) 
% diff内依次是函数、自变量、微分阶数,方程用==表示相等而不是赋值

2.初值问题

求初值问题

s = dsolve(\'Dy = y - 2*t / y\',\'y(0) =1\');

3.边界问题

求边界问题

s = dsolve(\'x*D2y - 3*Dy =x^2\',\'y(1)=0\',\'y(5) = 0\',\'x\');

4.高阶方程

求解方程

s=dsolve(\'D2y =cos(2*x) - y\',\'y(0) =1\',\'Dy(0) = 0\',\'x\');
simplify(s);
(eqn,cond,‘IgnoreAnalyticConstraints’,false) %设置不化简结果

5.方程组问题

求解方程组

[f,g]= dsolve(\'Df = f + g\',\'Dg = -f + g\',\'f(0)=1\',\'g(0) = 2\',\'x\'); 

一些例子

 dsolve(\'D2y+4*Dy+29*y = 0\',\'y(0) = 0\',\'Dy(0)= 15 \',\'x\')
 
ans =
 
3*sin(5*x)*exp(-2*x)

[x y z] =   dsolve(\'Dx = 2*x-3*y+3*z\',\'Dy = 4*x-5*y+3*z\',\'Dz = 4*x-4*y+2*z\')
 
x = 
C7*exp(2*t) + C8*exp(-t)
y =
C7*exp(2*t) + C8*exp(-t) + C9*exp(-2*t) 
z =
C7*exp(2*t) + C9*exp(-2*t)
%可以对其进行简化操作
 x = simplify(x) 
x = C7*exp(2*t) + C8*exp(-t)
 y = simplify(y) 
y =exp(-2*t)*(C9 + C8*exp(t) + C7*exp(4*t))

数值解

%龙格库塔法(Runge-Kutta法)
xfun=@(t,x)0.3.*x.*(1-x/8); %定义赋值函数r=0.3,k=8
[tout,xout]=ode45(fun,[0,40],0.1)  %方程数值解,四五阶RK法
[tout,xout]=ode23(xfun,[t0,tfinal],x0) %二三阶RK法
%%
ode系列数值求解形如  /  =   ( , )的微分方程组, 并绘图。
xfun: 输入参数,函数必须恰有t,x两个变量,用函数文件定义的fun.m则用@fun或‘fun’调用。
t0:输入参数,t的初始值。
tfinal:输入参数,t的终值。
x0:输入参数,x的初始值。
tout: 离散的自变量值, xout: 离散的函数值。
%%

同时也有一些其他的求解语句和输出语句

%%
其他的求解语句
ode45        ode113        ode15s 
ode23s       ode23t        ode23tb 
其他的输出语句
odeplot     odeprint        
odephas2    odephas3
%%

一个例子

求的数值解

首先对该方程进行换元然后建立m文件

function fyy=rhf(t,x)
    fyy=[y(1).*(1-y(2).^2)+y(2);y(1)];
end

最后计算数值解

y0=[0.25,0]’;
[t,y]=ode23(‘rhf’,[0,0.25],y0);
plot(t,y)

一些例子

%vdp1000.m

function dy = vdp1000(t,y)
dy = zeros(2,1);
dy(1) = y(2);
dy(2) = 1000*(1-y^2)*y2-y1;
end
%命令行输入
[T,Y] = ode15s(\'vdp1000\',[0 3000],[2 0]);%第一个参数是文件名,第二个参数是初始时间和终止时间第三个参数是y1和y2的初值
plot(T,Y(:,1),\'-\');
%结果是T时间 

plot(T,Y(:,1),\'-k\');,画Y数组中的第一列数随着T的变化曲线,‘-k’表示颜色黑色实线,

%定义函数
function dy=eq1(x,y)
    dy=zeros(2,1);
    dy(1)=y(2);
    dy(2)=1/5*sqrt(1+y(1)^2)/(1-x);
end

调用
x0=0;
xf=0.9999;
[x,y]=ode15s(\'eq1\',[x0 xf],[0 0]);
plot(x,y(:,1),\'-\') 
hold on
y=0:0.01:2;
plot(1,y,\'*\')

微分方程模型

1.种群增长Logistic模型

  • N(t)表示在时刻 t时刻种群数量
  • r 表示种群的内禀增长率,即在没有资源限制下的种群增长率
  • K表示环境载量,反映资源环境对种群增长的制约作用

2.生物种群竞争模型

  • N1(t)N2(t) 分别表示在时刻 甲、乙两个种群数量。
  • a11 表示种群甲自身的被抑制的情况
  • a12 表示种群乙对种群甲的竞争力

参考: https://zhuanlan.zhihu.com/p/162296418


鲜花

握手

雷人

路过

鸡蛋
该文章已有0人参与评论

请发表评论

全部评论

专题导读
上一篇:
delphi向另一程序窗口某处发送鼠标事件所用的工具发布时间:2022-07-18
下一篇:
Python for delphi教程发布时间:2022-07-18
热门推荐
阅读排行榜

扫描微信二维码

查看手机版网站

随时了解更新最新资讯

139-2527-9053

在线客服(服务时间 9:00~18:00)

在线QQ客服
地址:深圳市南山区西丽大学城创智工业园
电邮:jeky_zhao#qq.com
移动电话:139-2527-9053

Powered by 互联科技 X3.4© 2001-2213 极客世界.|Sitemap