- As I see it, the problem is at t = 0, y(1,2) = 0
ODE solver - division by zero at time boundaries
10 views (last 30 days)
Show older comments
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 Comments
Accepted 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?
0 Comments
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!