- As I see it, the problem is at t = 0, y(1,2) = 0
ODE solver - division by zero at time boundaries
9 vues (au cours des 30 derniers jours)
Afficher commentaires plus anciens
Bastian Andersen
le 12 Déc 2019
Réponse apportée : darova
le 16 Déc 2019
Hi
I have a problem when solving my system of ODEs. I simplified it to what I find essential for the problem.
As I see it, the problem is at t = 0, y(1,2) = 0 , whereas dy(2:end,2) = NaN AND at t = p.t_empty, y(1,1) = 0, whereas dy(2:end,1) = NaN
This is however the time interval I am interested in the solution for. It is essentially a system of two tanks. A reaction occurs in tank 1 and the matter is transferred to tank 2 meanwhile.
Do you have any suggestions to how I can solve this?
p.Q = 40;
p.Q_R = 150;
p.V = 20;
p.n = 101;
p.t_empty = p.V/p.Q;
y0 = zeros(p.n*2+2,1);
y0(1:2) = [p.V, 1];
options = odeset('RelTol',1e-12,'AbsTol',1e-12);
[t,y] = ode45('ODE_tank',[0,p.t_empty],y0,options,p);
figure(1)
plot((0:p.n-1),y(end,2:p.n+1))
figure(2)
plot((0:p.n-1),y(end,p.n+3:end))
function dy = ODE_tank(t,y,options,p)
y = reshape(y,[],2);
dy = zeros(size(y));
dy(1,1) = -p.Q;
dy(2,1) = -p.Q_R/y(1,1)*y(2,1);
dy(3:p.n,1) = p.Q_R/y(1,1)*(y(2:p.n-1,1) - y(3:p.n,1));
dy(p.n+1,1) = p.Q_R/y(1,1)*y(p.n,1);
dy(1,2) = p.Q;
dy(2:end,2) = p.Q/y(1,2)*(y(2:end,1)-y(2:end,2));
dy = reshape(dy,[],1);
end
2 commentaires
Réponse acceptée
darova
le 16 Déc 2019
- As I see it, the problem is at t = 0, y(1,2) = 0
Can you replace 0 with 1e-3?
0 commentaires
Plus de réponses (0)
Voir également
Catégories
En savoir plus sur Ordinary Differential Equations dans Help Center et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!