matlab微分方程转化状态方程_matlab求微分方程的通解步骤

matlab微分方程转化状态方程_matlab求微分方程的通解步骤[转]http://blog.sina.com.cn/s/blog_46e9b2010100tsqv.html用matlab时间也不短了,可是一直没有接触过微分方程。这次看看书,学习学习,记点儿笔记。

大家好,又见面了,我是你们的朋友全栈君。如果您正在找激活码,请点击查看最新教程,关注关注公众号 “全栈程序员社区” 获取激活教程,可能之前旧版本教程已经失效.最新Idea2022.1教程亲测有效,一键激活。

Jetbrains全系列IDE使用 1年只要46元 售后保障 童叟无欺

[转] http://blog.sina.com.cn/s/blog_46e9b2010100tsqv.html

 

用matlab时间也不短了,可是一直没有接触过微分方程。这次看看书,学习学习,记点儿笔记。

1.可以解析求解的微分方程。
dsolve()
调用格式为:

y=dsolve(f1,f2,…,fmO;

y=dsolve(f1,f2,…,fm,’x’);

如下面的例子,求解了微分方程
eq1.jpg
syms t;
u=exp(-5*t)*cos(2*t-1)+5;
uu=5*diff(u,t,2)+4*diff(u,t)+2*u;
syms t y;
y=dsolve([‘D4y+10*D3y+35*D2y+50*Dy+24*y=87*exp(-5*t)*cos(2*t-1)+92*exp(-5*t)*sin(2*t-1)+10’])
yc=latex(y)

将yc的内容copy到latex中编译,得到结果。

关于Matlab的微分方程,直到今天才更新第2篇,实在是很惭愧的事——因为原因都在于太懒惰,而不是其他的什么。

在上一篇中,我们使用dsolve可以解决一部分能够解析求解的微分方程、微分方程组,但是对于大多数微分方程(组)而言不能得到解析解,这时数值求解也就是没有办法的办法了,好在数值解也有很多的用处。

数值分析方法中讲解了一些Eular法、 Runge-Kutta 法等一些方法,在matlab中内置的ode求解器可以实现不同求解方法的相同格式的调用,而不必太关心matlab究竟是用什么算法完成的。

这一回我们来说明ode45求解器的使用方法。

1.ode45求解的上手例子:

求解方程组

Dx=y+x(1-x^2-y^2);

Dy=-x+y*(1-x^2-y^2)

初值x=0.1;y=0.2;

 

先说明一下最常用的ode45调用方式,和相应的函数文件定义格式。

[t,x]=ode45(odefun,tspan,x0);

其中,Fun就是导函数,tspan为求解的时间区间(或时间序列,如果采用时间序列,则必须单调),x0为初值。

这时,函数文件可以采用如下方式定义

function dx=odefun(t,x)

对于上面的小例子,可以用如下的程序求解。

 

 

function jixianhuan
clear;clc
x0=[0.1;0.2];
[t,x]=ode45(@jxhdot,[0,100],x0);
plot(x(:,1),x(:,2))

 

function dx=jxhdot(t,x)
dx=[
  x(2)+x(1).*(1-x(1).^2-x(2).^2);
  -x(1)+x(2).*(1-x(1).^2-x(2).^2)
];

 

 

 

初识Matlab微分方程(2)

 2.终值问题

tspan可以是递增序列,也可以为递减序列,若为递减则可求解终值问题。

初识Matlab微分方程(2)初识Matlab微分方程(2)

[t,x]=ode45(@zhongzhiode,[3,0],[1;0;2]);plot(t,x)

 

function dx=zhongzhiode(t,x)
dx=[2*x(2)^2-2;
-x(1)+2*x(2)*x(3)-1;
-2*x(2)+2*x(3)^2-4];

结果如下

 

初识Matlab微分方程(2)

3.odeset

 

options = odeset(‘name1′,value1,’name2’,value2,…)

 

[t,x]=solver(@fun,tspan,x0,options)

通过odeset设置options

第一,通过求解选项的设置可以改善求解精度,使得原本可能不收敛的问题收敛。

options=odeset(‘RelTol’,1e-10);

第二,求解形如M(t,x)x’=f(t,x)的方程。

例如,方程

x’=-0.2x+yz+0.3xy

y’=2xy-5yz-2y^2

x+y+z-2=0

可以变形为

[1  0  0][x’]    [-0.2x+yz+0.3xy]

[0  1  0][y’] = [2xy-5yz-2y^2   ]

[0  0  0][z’]    [x+y+z-2           ]

这样就可以用如下的代码求解该方程

function mydae
M=[1 0 0;0 1 0;0 0 0];
options=odeset(‘Mass’,M);
x0=[1.6,0.3,0.1];
[t,x]=ode15s(@daedot,[0,1.5],x0,options);plot(t,x)

function dx=daedot(t,x)
dx=[
    -0.2*x(1)+x(2)*x(3)+0.3*x(1)*x(2);
    2*x(1)*x(2)-5*x(2)*x(3)-2*x(2)*x(2);
    x(1)+x(2)+x(3)-2];
初识Matlab微分方程(3)

 

4.带附加参数的ode45

有时我们需要研究微分方程组中的参数对于解的影响,这时采用带有参数的ode45求解会使求解、配合循环使用,可以使得求解的过程更加简捷。

使用方法:只需将附加参数放在options的后面就可以传递给odefun了。

看下面的例子。

function Rossler

clear;clc
a=[0.2,0.2];
b=[0.2,0.5];
c=[5.7,10];

x0=[0 0 0];
for jj=1:2
    [t,x]=ode45(@myRossler,[0,100],x0,[],a(jj),b(jj),c(jj));
    figure;plot3(x(:,1),x(:,2),x(:,3));grid on;
end

function dx=myRossler(t,x,a,b,c)

dx=[
    -x(2)-x(3);
    x(1)+a*x(2);
    b+(x(1)-c)*x(3)];
初识Matlab微分方程(3)初识Matlab微分方程(3)

5. 刚性方程的求解

刚性方程就是指各个自变量的变化率差异很大,会造成通常的求解方法失效。

这是matlab中自带的一个例子,使用ode15s求解,如果用ode45求解就会出现错误。

初识matlab微分方程(4)

function myode15study

[t,Y] = ode15s(@vdp1000,[0 3000],[2 0]);
plot(T,Y(:,1),’-o’)

figure;plot(Y(:,1),Y(:,2))

function dy = vdp1000(t,y)
dy = zeros(2,1);   

dy(1) = y(2);
dy(2) = 1000*(1 – y(1)^2)*y(2) – y(1);

初识matlab微分方程(4)

初识matlab微分方程(4)

 

6.高阶微分方程的求解

通常的方法是进行变量替换,将原方程降阶,转换成更多变量的一阶方程组进行求解。

在这个例子里我们求解一个动力学系统里最常见的一个运动方程

初识matlab微分方程(4),其中f=sin(t)

初识matlab微分方程(4)

 

function myhighoder
clear;clc
x0=zeros(6,1);
[t,x]=ode45(@myhigh,[0,100],x0);
plot(t,x(:,1))

function dx=myhigh(t,x)
f=[sin(t);0;0];;
M=eye(3);
C=eye(3)*0.1;
K=eye(3)-0.5*diag(ones(2,1),1)-0.5*diag(ones(2,1),-1);
dx=[x(4:6);inv(M)*(f-C*x(4:6)-K*x(1:3))];

 

7.延迟微分方程

matlab提供了dde23求解非中性微分方程。dde23的调用格式如下:

sol = dde23(ddefun,lags,history,tspan)

lags是延迟量,比如方程中包含y1(t-0.2)和y2(t-0.3)则可以使用lags=[0.2,0.3]。

这里的ddefun必须采用如下的定义方式:

dydt = ddefun(t,y,Z)

其中的Z(:,1)就是y(t-lags(1)),Z(:,2)就是y(t-lags(2))…

下面是个使用dde23求解延迟微分方程的例子。

function mydde23study
%   The differential equations
%
%        y’_1(t) = y_1(t-1) 
%        y’_2(t) = y_1(t-1)+y_2(t-0.2)
%        y’_3(t) = y_2(t)
%
%   are solved on [0, 5] with history y_1(t) = 1, y_2(t) = 1, y_3(t) = 1 for
%   t <= 0.
clear;clc
lags=[1,0.2];
history=[1;1;1];
tspan=[0,5];
sol = dde23(@myddefun,lags,history,tspan)
plot(sol.x,sol.y)

function dy = myddefun(t,y,Z)
dy=[
    Z(1,1);
    Z(1)+Z(2,2);
    y(2)    ];

 

初识matlab微分方程(5)

 

8.ode15i求解隐式微分方程

[T,Y] = ode15i(odefun,tspan,y0,yp0)

yp0为y’的初值。

odefun的格式如下  dy = odefun(t,y,yp),yp表示y’,而方程中应该使得f(t,y,y’)=0

 

function myodeIMP
%   The problem is
%
%         y(1)’ = -0.04*y(1) + 1e4*y(2)*y(3)
%         y(2)’ =  0.04*y(1) – 1e4*y(2)*y(3) – 3e7*y(2)^2
%         y(3)’ =  3e7*y(2)^2
%
%   It is to be solved with initial conditions y(1) = 1, y(2) = 0, y(3) = 0
%   to steady state. 
clear;clc
y0=[1;0;0];
fixed_y0=[1;1;1];
yp0=[0 0 0];
fixed_yp0=[];

[y0mod,yp0mod]=decic(@myodefunimp,0,y0,fixed_y0,yp0,fixed_yp0);
tspan=[0, logspace(-6,6)];
[t,y] = ode15i(@myodefunimp,tspan,y0mod,yp0mod);
y(:,2)=1e4*y(:,2);
semilogx(t,y)
function res=myodefunimp(t,y,yp)
res=[
    -yp(1)-0.04*y(1)+1e4*y(2)*y(3);
    -yp(2)+0.04*y(1)-1e4*y(2)*y(3)-3e7*y(2)^2;
    -yp(3)+3e7*y(2)^2;
    ];
初识matlab微分方程(5)

这次要接触一个新的求解ode的方法,就是使用simulink的积分器求解。

1.还是做我们研究过的一个例子(在初识matlab微分方程(2)中采用的)。

Dx=y+x(1-x^2-y^2);

Dy=-x+y*(1-x^2-y^2)

初值x=0.1;y=0.2;

积分器中设置初始条件;f(u)中指定Dx,Dy的计算公式。初识Matlab微分方程(6)

运行这个仿真,scope中可以看到两个变量的时程如下:

初识Matlab微分方程(6)

在WorkSpace里可以得到tout和yout,执行plot(yout(:,1),yout(:,1))得到与ode45求解相似的结果如下

初识Matlab微分方程(6)

 2.这部分解决一个使用ode求解器dde23没法求解的一类延迟微分方程(中性微分方程)。

形如x'(t)=f(x'(t-t1),x(t),x(t-t2),x(t-t3))这类方程。dde23是无法求解的,但是可以借助simulink仿真求解。

看下面的这个例子。

x'(t)=A1*x(t-t1)+A2*x'(t-t2)+B*u(t)

t1=0.15;t2=0.5

A1=[-12     3   -3]      A2=[0.02    0     0]    B=[0]

   [106  -116   62]         [0    0.03     0]      [1]

   [207  -207  113]         [0       0  0.04]      [2]   

在continuous里找到transport Delay,就可以实现对于信号的延迟,因此可以建立如下仿真模型

初识Matlab微分方程(6)

从而在scope中可以得到如下仿真结果

初识Matlab微分方程(6)

 

 

OK~初识微分方程到了这里我想应该可以做个终结,因为我想作为零基础的材料来看,到这里也就可以了。以后还可能再有微分方程的内容,还请感兴趣的朋友多捧场吧。

最后,大力推荐一本书薛定宇老师的《高等应用数学问题的Matlab求解》,确实很经典。学习Matlab的时间也不算短了,可是每次翻看这本书总是能让我有温故而知新的感觉,是我目前见过的最好的Matlab书。强烈推荐!(对于从来没有接触过matlab的人来说或许有点儿难,但是如果你以后要用matlab的话买一本绝对不会后悔的。)

版权声明:本文内容由互联网用户自发贡献,该文观点仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 举报,一经查实,本站将立刻删除。

发布者:全栈程序员-用户IM,转载请注明出处:https://javaforall.cn/167738.html原文链接:https://javaforall.cn

【正版授权,激活自己账号】: Jetbrains全家桶Ide使用,1年售后保障,每天仅需1毛

【官方授权 正版激活】: 官方授权 正版激活 支持Jetbrains家族下所有IDE 使用个人JB账号...

(0)
blank

相关推荐

  • JavaWeb之HttpSession

    JavaWeb之HttpSessionHttpSession一、概述HttpSession是由JavaWeb提供的,用来会话跟踪的类。session是服务器端对象,保存在服务器端!!!HttpSession是Servlet三大域对象之一,所以它也有setAttribute()、getAttribute()、removeAttribute()方法。HttpSession底层依赖Cookie,或是URL重写!二、HttpSe…

  • porm文件中的build标签

    porm文件中的build标签分类    (1)全局配置(projectbuild)         针对整个项目的所有情况都有效    (2)配置(profilebuild)         针对不同的profile配置[html] viewplain copyproject xm

  • mybatis异常invalid comparison: java.util.Date and java.lang.String

    mybatis异常invalid comparison: java.util.Date and java.lang.String开发中改动mapper文件后需要重新编译发布,由于工程比较大非常耗时,所以为方便快速测试干脆写了一个小java工程.工程中用到的dao,mapper和实体类都是从工程中拷出来的,数据库也是同一个.但是遇到一个比较奇怪的问题实体类中有一个属性privateDatecreateTime;对应该属性数据库中定义的是create_timedatetime

  • 决不随波逐流,记住15件事让你达成目标

    决不随波逐流,记住15件事让你达成目标

    2021年11月29日
  • Android程序员接私活完整攻略「建议收藏」

    Android程序员接私活完整攻略「建议收藏」接私活对程序员这个圈子来说是一个既公开又隐私的话题,不说全部,应该大多数程序员都有过想要接私活的想法,当然,也有部分得道成仙的不主张接私活。但是很少有人在公开场合讨论私活的问题,似乎都在避嫌。就跟有人下班后跑滴滴一样,程序员私有时间接点活挣点钱不也很正常么,不过不要在上班时间就行,就跟你上班期间出去跑滴滴一样。当你竭尽全力想要去接私活的时候一定做过这样的事,百度搜索“程序员如何接私活”或者…

  • Kong集群(hybrid混合)部署模式

    Kong集群(hybrid混合)部署模式kong集群部署与promethues指标采集

发表回复

您的电子邮箱地址不会被公开。

关注全栈程序员社区公众号