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
采纳的回答
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 个评论
@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
% 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
Torsten
2023-8-28
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
2023-8-28
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')
更多回答（0 个）
类别
在 帮助中心 和 File Exchange 中查找有关 Mathematics 的更多信息
Translated by 
