# ODE solver - division by zero at time boundaries

13 views (last 30 days)
Bastian Andersen on 12 Dec 2019
Answered: darova on 16 Dec 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 CommentsShowHide 1 older comment
Bastian Andersen on 16 Dec 2019
This works, thank you. How do I accept this as the answer?

darova on 16 Dec 2019
• As I see it, the problem is at t = 0, y(1,2) = 0
Can you replace 0 with 1e-3?

### Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by