隐式微分方程求解Matlab

微分方程求解Matlab详解

   在matlab求解微分方程有很多中，ode45是最常见的，当然，ode45求解有一定条件，就是需要把高阶微分方程转化为1阶微分方程。即使用降阶法。降阶法原理就是引入中间变量，然后实现降阶。例如：$$F(t,y,y^{{}^{\prime }},y^{{}^{\prime \prime }})=0$$   令$x_1=y$, $x_2=y^{{}^{\prime }}$,那么就可以将原来的二阶微分方程降阶为两个一阶微分方程。
$$\begin{array}{c}\frac{d{x}_1}{dt}=x_2\\ \frac{d{x}_2}{dt}=G(t,x_1,x_2)\end{array}$$
   高阶微分方程组同理！！！

例子

   这里以一个简单的例子演示ode45，以及包括之后使用的均为这个例子：
$$y^{{}^{\prime }}-2t=0$$
   并且指定时间区间 $[0,5]$和初始条件 $y_0=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}^{{}^{\prime \prime }3}+y^{{}^{\prime }}{y}^{{}^{\prime \prime }}+y=0$$
   那对于这种该如何求解呢？

方法1-使用solve()求解高阶表达式

   以上面的例子为例，将$F(t,y,y^{{}^{\prime }})=0$看做是$y^{{}^{\prime }}$的方程，利用solve()函数求解出$y^{{}^{\prime }}$的表达式。这种方法是省去了手动的降阶，而使用计算机代替我们计算，我将这个封装了一个函数，可以通用，源码和参考在这里；

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^{{}^{\prime }})=0$，输入是t,y，输出是$y^{{}^{\prime }}$,也就是说，我们要通过变量的低阶值和时间t求解变量的最高阶的值，因此，对于任意复杂的$F(t,y,y^{{}^{\prime }})=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')

源码gtihub

    GitHub传送门！！！