ODE45 running infinitely without solving

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 !

Answer 1

[Walter Roberson]

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

[Ajmal R S]

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]

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]

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]

we need the rest of your code to test with

[Ajmal R S]

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