MATLAB solve ODE on invariant manifold

Question

I have a system that looks like

dn/dt=f(n,v)
dh/dt=g(h,v)

I want to solve this equation on the manifold F(v,n,h)=0, a nonlinear function in v. I tried to use something like v=fzero(@(x) F(x,n,h),0) to solve for the value of v on the manifold at each time step. But this is incredibly slow and ode15s (my system is a relaxation oscillator) fails to meet integration tolerance. How do I find solution to the ODE on the manifold defined by F(v,n,h)=0?

You should employ a DAE solver, as your system is differential-algebraic. — Lutz Lehmann, Oct 03 '16 at 23:54

Badoe · Accepted Answer · 2016-10-08T16:19:56.483

I find @LutzL's comment very helpful. It is possible to set up a DAE solver using ode15s. Example: The "Solve Robertson Problem as Semi-Explicit Differential Algebraic Equations (DAEs)" section in https://www.mathworks.com/help/matlab/ref/ode15s.html

In my case, I would setup a matrix:

M=[zeros(1,3);0,1,0;0,0,1];
options = odeset('Mass',M,'RelTol',1e-5,'AbsTol',1e-6,'MaxStep',0.01);
y0=[v0,n0,h0];
[T,Y]=ode15s(@slow,[0 50],y0,options);

And slow is a function defined as:

function dy = slow(t,y)
    v=y(1); n=y(2); h=y(3);
    dy=zeros(3,1);
    dy(1)=F(v,n,h);
    dy(2)=f(n,v);
    dy(3)=g(h,v);
end

score 0 · Answer 2 · edited Feb 09 '17 at 15:11

0

One of the possible ways to solve this problem is to differentiate the F(v,n,h)=0 equation:

$\frac{\partial F}{\partial v}\frac{dv}{dt}+\frac{\partial F}{\partial n}\frac{dn}{dt}+\frac{\partial F}{\partial h}\frac{dh}{dt}=0$

Now we can obtain the ODE system

$\frac{dn}{dt}=f(n,v),\, \frac{dh}{dt}=g(h,v),\,\frac{dv}{dt}= -\left.\left( \frac{\partial F}{\partial n}\frac{dn}{dt}+\frac{\partial F}{\partial h}\frac{dh}{dt} \right )\right/ \left( \frac{\partial F}{\partial v} \right )$

or

$\frac{dn}{dt}=f(n,v),\, \frac{dh}{dt}=g(h,v),\,\frac{dv}{dt}= -\left.\left( \frac{\partial F}{\partial n}f(n,v)+\frac{\partial F}{\partial h}g(h,v) \right )\right/ \left( \frac{\partial F}{\partial v} \right )$

which can be solved in the usual way.

edited Feb 09 '17 at 15:11

Community

1
1

answered Oct 04 '16 at 03:05

AVK

2,130
3
15
25

I have tried this way. However, the dF/dv term becomes very close to 0 near a fold structure in the system. Hence the solution precision is not very good. – Badoe Oct 04 '16 at 14:02

MATLAB solve ODE on invariant manifold

2 Answers2