解微分方程方法

MATLAB解微分方程(2011-07-15 17:35:25)

转载▼

分类：matlab学习标签：

教育

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

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

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

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

function dx=odefun(t,x)

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

2.终值问题

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

[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];

结果如下

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 1][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];

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)];

5. 刚性方程的求解

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

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

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);

6.高阶微分方程的求解

通常的方法是进行变量替换，将原方程降阶，转换成更多变量的一阶方程组进行求解。在这个例子里我们求解一个动力学系统里最常见的一个运动方程

,其中f=sin(t)

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

%