ODE45 running infinitely without solving

2020 5 4
1 个回答
回答已采纳
更新时间：2020 5 5
2 次查看（30 天）

0 个投票

I tried to code an ode45 function to find the concentration wrt to thickness of membrane but the ode45 function is taking too long without actually solving the code
ts = 100e-06; %in m
ta = 100e-09; %in m
z0 = ta + ts ; %z = ta+ts
zend = ta + ts + del_d ; %z = ta+ts+delta_d
[zdsol,cdsol]=ode45(@(z,c) decp(z,c), [z0,zend],c0);
And this is the function
function dc = decp(z,c)
x1_dsd = -6.92e-12;
x2_dsd = -9.95e-11;
x3_dsd = 1.51e-9;
dsd_d= x1_dsd*c^2 + x2_dsd*c + x3_dsd;
dc = ((0.02+6*c)/dsd_d);
end
With Js and Jw being constants
Please help !

请先登录，再回答此问题。

 采纳的回答

Walter Roberson
Walter Roberson 2020-5-5
编辑：Walter Roberson 2020-5-5

0 个投票

try ode23s. if your system is stiff then ode45 can take a long long time.
Also I worry that your c0 is associated with time 0; it needs to be associated with time z0 instead.
Your denominator dsd_c has zeros near time 1e-11 which could be a problem

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Ajmal R S
Ajmal R S 2020-5-5
Thank you,
I corrected the c0 problem and changed it to ode23s, but now i am getting negative values beyond the first few solutions.
Is it because of my dsd_c ?
Walter Roberson
Walter Roberson 2020-5-5
we need enough to be able to repeat the problem.
But your slope is dominated by negative times c squared. that is going to be negative unless your initial conditions are large enough. But if it is large then at some point in decreasing it would pass through a zero in the denominator. There might perhaps be an area near 0 where the equation causes oscillating between negative and positive
Ajmal R S
Ajmal R S 2020-5-5
I changed the coefficient of c squared to positive, but still I am getting negative values, the rest of code mostly deals with finding the value of del_d. The value of del_d is 7.637240610884750e-05. The integration is with z0= del_d and zend= 0 with c0=1 corresponding to value of c at z=z0.
Walter Roberson
Walter Roberson 2020-5-5
we need the rest of your code to test with
Ajmal R S
Ajmal R S 2020-5-5
Here is the full code.
clc;
clear all;
A = 0.5;
B = 0.4;
S = 200e-06;
Jw= 6;
Js= 0.02;
c0 = 1;
v = 20; %cross flow velocity feed and draw in cm/s
x1_rho = 0.04;
x2_rho = 1;
rho_d = x1_rho*c0 + x2_rho; %in cgs
x1_mu = 0.02;
x2_mu = 0.06;
x3_mu = 0.89;
mu_d = x1_mu*c0^2 + x2_mu*c0 + x3_mu; %in cgs
L = 8.55; %channel length in cm
dh = 0.42; %hydrauic diameter in cm made it 0.2 to meet conditions in ref.
Re = (rho_d*1000*v*10^-2*dh*10^-2)/(mu_d*10^-3);
x1_dsd = 6.92e-12;
x2_dsd = -9.95e-11;
x3_dsd = 1.51e-9;
dsd_d= x1_dsd*c0^2 + x2_dsd*c0 + x3_dsd; %given in section 3.2 in m2/s
Sc = (mu_d*10^-3)/(rho_d*1000*dsd_d);
Sh = 1.85*(Re*Sc*(dh/L))^0.33;
k_d = (Sh*dsd_d)/(dh*10^-2); %in m/s
del_d = dsd_d/k_d; %in m
ts = 100e-06; %in m
ta = 100e-09; %in m
zend = ta + ts ; %z = ta+ts
z0 = ta + ts + del_d ; %z = ta+ts+delta_d
[zdsol,cdsol]=ode23s(@(z,c) decp(z,c), [z0,zend],c0);
And the function
function dc = decp(z,c)
x1_dsd = 6.92e-12;
x2_dsd = -9.95e-11;
x3_dsd = 1.51e-9;
dsd_d= x1_dsd*c^2 + x2_dsd*c + x3_dsd;
dc = ((0.02+6*c)/dsd_d);
end

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息

产品

版本

R2019b

标签

Community Treasure Hunt

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

Start Hunting!

Translated by