微分方程求解Matlab详解
在matlab求解微分方程有很多中,ode45是最常见的,当然,ode45求解有一定条件,就是需要把高阶微分方程转化为1阶微分方程。即使用降阶法。降阶法原理就是引入中间变量,然后实现降阶。例如:
F
(
t
,
y
,
y
′
,
y
′
′
)
=
0
F(t,y,y^{'},y^{''})=0
F(t,y,y′,y′′)=0 令
x
1
=
y
x_1=y
x1=y,
x
2
=
y
′
x_2=y^{'}
x2=y′,那么就可以将原来的二阶微分方程降阶为两个一阶微分方程。
d
x
1
d
t
=
x
2
d
x
2
d
t
=
G
(
t
,
x
1
,
x
2
)
\begin{matrix} \frac{dx_1}{dt}=x_2 \\ \frac{dx_2}{dt}=G(t,x_1,x_2) \end{matrix}
dtdx1=x2dtdx2=G(t,x1,x2)
高阶微分方程组同理!!!
例子
这里以一个简单的例子演示ode45,以及包括之后使用的均为这个例子:
y
′
−
2
t
=
0
y^{'}-2t = 0
y′−2t=0
并且指定时间区间
[
0
,
5
]
[0,5]
[0,5]和初始条件
y
0
=
0
y_{0} = 0
y0=0。
ODE45求解
使用匿名函数的方式
clc;clear all;
tspan = [0 5];
y0 = 0;
figure(1)
[t,y] = ode45(@(t,y) 2*t, tspan, y0);
plot(t,y,'-o')
使用.m文件的方式
function dy = Fcn(t,y)
%FCN
dy=2*t;
end
clc;clear all;
tspan = [0 5];
y0 = 0;
figure(2)
[t1,y2] = ode45(@Fcn, tspan, y0);
plot(t1,y2,'-o')
隐式微分方程求解
ODE45使用最困难的一步就是微分方程的降阶,对于一个含有最高阶耦合项的微分方程来说,根本无法手动将其化简为一阶微分方程,例如:
t
y
′
′
3
+
y
′
y
′
′
+
y
=
0
ty^{''3}+y^{'}y^{''}+y = 0
ty′′3+y′y′′+y=0
那对于这种该如何求解呢?
方法1-使用solve()求解高阶表达式
以上面的例子为例,将 F ( t , y , y ′ ) = 0 F(t,y,y^{'})=0 F(t,y,y′)=0看做是 y ′ y^{'} y′的方程,利用solve()函数求解出 y ′ y^{'} y′的表达式。这种方法是省去了手动的降阶,而使用计算机代替我们计算,我将这个封装了一个函数,可以通用,源码和参考在这里;
clc;clear all;
syms dy t y
eq = dy-2*t;
dy = solve(eq,dy)
matlabFunction(dy,'File','Fcn','Vars',[t,y]);
tspan = [0 5];
y0 = 0;
figure(1)
[t,y] = ode45(@Fcn, tspan, y0);
plot(t,y,'-o')
方法2-使用求解高阶变量的数值解,再使用ode45
然而,对于上述的方法,不是所有的微分方程都能算出最高阶的解析解。那么我们就么有办法了吗?观察ode45所传入的函数句柄,例如 F ( t , y , y ′ ) = 0 F(t,y,y^{'})=0 F(t,y,y′)=0,输入是t,y,输出是 y ′ y^{'} y′,也就是说,我们要通过变量的低阶值和时间t求解变量的最高阶的值,因此,对于任意复杂的 F ( t , y , y ′ ) = 0 F(t,y,y^{'})=0 F(t,y,y′)=0或者是高阶微分方程组,我们在函数句柄,利用fsolve求解出最高阶的值,将其返回即可。
function dy = Fcn(t,y)
fun = @(dy)dy-2*t;
options=optimset('display','off');
dy = fsolve(fun,0,options) ;
clc;clear all;
tspan = [0 5];
y0 = 0;
figure(2)
[t1,y2] = ode45(@Fcn, tspan, y0);
plot(t1,y2,'-o')
方法3-使用ode15i
matlab在后续版本提供了ode15i,其求解思路类似上面方法2。具体使用:
function fcn = Fcn(t,y,dy)
%WEISSINGER Evaluate the residual of the Weissinger implicit ODE
%
% See also ODE15I.
% Jacek Kierzenka and Lawrence F. Shampine
% Copyright 1984-2014 The MathWorks, Inc.
fcn = dy - t.*2.0;
clc;clear all;
t0 = 0;
tspan = [t0 5];
y0 = 0;
yp0 =1;
% 计算一致的初始条件
% ode15i 求解器需要一致的初始条件,即提供给求解器的初始条件必须满足
% 这个就是需要满足微分方程
[y0,dy0] = decic(@Fcn,t0,y0,1,yp0,0);
[t,y] = ode15i(@Fcn, tspan, y0,dy0);
figure(1)
plot(t,y,'-o')