Ordinary differential equations with two unknows (ode45)
2 views (last 30 days)
Show older comments
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
0 Comments
Accepted Answer
Sam Chak
on 28 Aug 2023
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
% 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
4 Comments
Sam Chak
on 29 Aug 2023
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')
More Answers (0)
See Also
Categories
Find more on Ordinary Differential Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!