フィルターのクリア

ODE solver - division by zero at time boundaries

9 ビュー (過去 30 日間)
Bastian Andersen
Bastian Andersen 2019 年 12 月 12 日
回答済み: darova 2019 年 12 月 16 日
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 件のコメント
darova
darova 2019 年 12 月 13 日
  • As I see it, the problem is at t = 0, y(1,2) = 0
Can you replace 0 with 1e-3?
Bastian Andersen
Bastian Andersen 2019 年 12 月 16 日
This works, thank you. How do I accept this as the answer?

サインインしてコメントする。

採用された回答

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

その他の回答 (0 件)

カテゴリ

Help Center および File ExchangeOrdinary Differential Equations についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by