Ordinary differential equations with two unknows (ode45)

Hi, thank you very much for any help!

I am trying to solve two equations by converting them into two ordinary differential equations, but ODE45, ODE23, et al. cannot solve them well.

The two unknows vv and ff are a function of time t. The code is listed as below:

p.a=0.0062; p.b=0.0128; p.L=17.5e-6; p.Vo=1e-6; p.Vl=5e-7; p.sigma=5e7; p.k=0.95*p.sigma*(p.b-p.a)/p.L;

odefun=@(t,y) ode_isothermal(p,t,y);
yo=zeros(2,1);
yo(1) = p.L/p.Vo;   % initial state variable
yo(2) = p.Vo;       % initial velocity

tic  % four/fifth order Runge-Kutta solver
options=odeset('Refine',1,'RelTol',1e-6,'InitialStep',1e-5,'MaxStep',1e5); 
[t,y]=ode45(odefun,[0 2000],yo,options);
toc

figure(1);  % figure
plot(t,log(y(:,2)/p.Vo));
ylabel('ln(V/V_{0})');xlable('Time');
box on; set(gca,'xlim',[0 2000])

%function
function dydt=ode_isothermal(p,~,y) 
dydt=zeros(2,1);
ff=y(1);
vv=y(2);

dydt(1)=1-ff*vv/p.L;             % function 1 for ff
kv=p.k*(p.Vl-vv);
bb=p.b*(dydt(1)/ff)*(p.Vl/vv).^(1/3);
aa=p.b/3/vv*(p.Vl/vv).^(1/3).*log(ff/35);
dydt(2)=(kv-bb)/(p.a/vv+0.1-aa);   %function 2 for vv
end

Answer 1

[Sam Chak]

Hi @Mei Cheng

If ode15s() solver is used, the simulation runs. You should be able to plot the results.

p.a     = 0.0062; 
p.b     = 0.0128; 
p.L     = 17.5e-6; 
p.Vo    = 1e-6; 
p.Vl    = 5e-7; 
p.sigma = 5e7; 
p.k     = 0.95*p.sigma*(p.b - p.a)/p.L;

odefun = @(t, y) ode_isothermal(p,t,y);
y0    = zeros(2, 1);
y0(1) = p.L/p.Vo;   % initial state variable
y0(2) = p.Vo;       % initial velocity

tic  % four/fifth order Runge-Kutta solver
% options = odeset('Refine', 1, 'RelTol', 1e-6, 'InitialStep', 1e-5, 'MaxStep', 1e5); 
[t, y] = ode15s(odefun, [0 2000], y0);
toc

Elapsed time is 0.107397 seconds.

% test
subplot(211)
plot(t, y(:,1)), grid on, xlabel('Time'), title('ff')
subplot(212)
plot(t, y(:,2)), grid on, xlabel('Time'), title('vv')

% figure(1);  % figure
% plot(t, log(y(:,2)/p.Vo));
% xlabel('Time'), ylabel('ln(V/V_{0})')
% box on; set(gca,'xlim',[0 2000])

%function
function dydt = ode_isothermal(p,~,y) 
    dydt    = zeros(2, 1);
    ff      = y(1);
    vv      = y(2);
    
    dydt(1) = 1 - ff*vv/p.L;             % function 1 for ff
    kv      = p.k*(p.Vl-vv);
    bb      = p.b*(dydt(1)/ff)*(p.Vl/vv).^(1/3);
    aa      = p.b/3/vv*(p.Vl/vv).^(1/3).*log(ff/35);
    dydt(2) = (kv - bb)/(p.a/vv + 0.1 - aa);   %function 2 for vv
end

[Mei Cheng]

@Sam Chak Thank you very much for your answering. I apologied for one mistake, that is, the value of p.k is just corrected. If I use p.k = 5*(p.b - p.a)/p.L, it is okay; but if I use p.k = 0.95*(p.b - p.a)/p.L, there is some problem. I don't know how to fix it.

p.a     = 0.0062; 
p.b     = 0.0128; 
p.L     = 17.5e-6; 
p.Vo    = 1e-6; 
p.Vl    = 5e-7; 
p.sigma = 5e7; 
p.k     = 0.95*(p.b - p.a)/p.L;

odefun = @(t, y) ode_isothermal(p,t,y);
y0    = zeros(2, 1);
y0(1) = p.L/p.Vo;   % initial state variable
y0(2) = p.Vo;       % initial velocity

tic  % four/fifth order Runge-Kutta solver
% options = odeset('Refine', 1, 'RelTol', 1e-6, 'InitialStep', 1e-5, 'MaxStep', 1e5); 
[t, y] = ode15s(odefun, [0 2000], y0);
toc

Elapsed time is 0.230605 seconds.

% test
subplot(211)
plot(t, y(:,1)), grid on, xlabel('Time'), title('ff')

Warning: Imaginary parts of complex X and/or Y arguments ignored.

subplot(212)
plot(t, y(:,2)), grid on, xlabel('Time'), title('vv')

Warning: Imaginary parts of complex X and/or Y arguments ignored.

% figure(1);  % figure
% plot(t, log(y(:,2)/p.Vo));
% xlabel('Time'), ylabel('ln(V/V_{0})')
% box on; set(gca,'xlim',[0 2000])

%function
function dydt = ode_isothermal(p,~,y) 
    dydt    = zeros(2, 1);
    ff      = y(1);
    vv      = y(2);
    
    dydt(1) = 1 - ff*vv/p.L;             % function 1 for ff
    kv      = p.k*(p.Vl-vv);
    bb      = p.b*(dydt(1)/ff)*(p.Vl/vv).^(1/3);
    aa      = p.b/3/vv*(p.Vl/vv).^(1/3).*log(ff/35);
    dydt(2) = (kv - bb)/(p.a/vv + 0.1 - aa);   %function 2 for vv
end

[Torsten]

Most probably, y(1) or y(2) become negative which will make aa and thus dydt(2) complex-valued. This is caused either by the log(ff/35) or the (p.Vl/vv)^(1/3) term.

[Mei Cheng]

@TorstenThank you very much for your further reply. I am not sure about it, but i will try to simplify my equations. Thank you again!

[Sam Chak]

Hi @Mei Cheng

As explained by @Torsten, the logarithmic and cubic root terms are primarily the causes of the complex-valued results. You can observe the streamlines in this plot to visually understand how the trajectories of the states "ff" and "vv" evolve as they approach 0. Additionally, it's worth noting that the equations also include several instances of division-by-zero terms.

% parameters
p.a     = 0.0062; 
p.b     = 0.0128; 
p.L     = 17.5e-6; 
p.Vo    = 1e-6; 
p.Vl    = 5e-7; 
p.sigma = 5e7; 
p.k     = 0.95*(p.b - p.a)/p.L;

% meshing
f       = linspace(-1.2, 1.2, 15);
v       = linspace(-0.3, 0.3, 15);
[ff,vv] = meshgrid(f, v);

% streamlines
dy1     = 1 - ff.*vv/p.L;
kv      = p.k*(p.Vl - vv);
bb      = p.b*(dy1./ff).*(p.Vl./vv).^(1/3);
aa      = p.b/3./vv.*(p.Vl./vv).^(1/3).*log(ff/35);
dy2     = (kv - bb)./(p.a./vv + 0.1 - aa);
l       = streamslice(ff, vv, dy1, dy2);

% axis control
axis tight
axis([-1, 1, -0.25, 0.25])
xlabel('ff'), ylabel('vv')