0
$\begingroup$

Does non-dimensionalizing ODEs affect the stiffness of the equations? Can it improve the stability of numerical methods like ode45,ode113 in MATLAB?

I am trying to solve 2 eqns. which might involve very tiny step size to get stability like 1e-11. I want to improve the speed of the numerical method. Can non dimensionalizing help in this regard?

The ODEs (I.V: Rbar, N) are shown in link 1. The output should be close to that shown in 2

My working code for a test problem:

clc
clear
close all
tic

T  = 273 + 160 ;  
Xo = 0.06;

[SS,delGv,rc,Z,rstar,delGstar,I_nu] = nucleate2(Xo, T, 0);

t  = [1 1e7];
y0 = [rstar,1];
options = odeset('InitialStep',1e-3,'MaxStep',100);
[t1,y1] = ode23(@(t,y)NOCmodel2(t,y,T),t,y0,options);

figure(1)
set(gcf, 'Position',  [100, 100, 900, 800])
loglog(t1,y1(:,1),'-b');
xlabel ('time(s)')
ylabel ('$R_{mean}$','Interpreter','latex')
set(findall(gcf,'-property','FontSize'),'FontSize',20,'fontweight','bold')


figure(2)
set(gcf, 'Position',  [100, 100, 900, 800])
loglog(t1,y1(:,2),'-r');

xlabel ('time(s)')
ylabel ('$logN (/mm^3)$','Interpreter','latex')
set(findall(gcf,'-property','FontSize'),'FontSize',20,'fontweight','bold')

figure(3)
set(gcf, 'Position',  [100, 100, 900, 800])
semilogx(t1,y1(:,2),'-r');
xlabel ('time(s)')
ylabel ('$N (/mm^3)$','Interpreter','latex')
set(findall(gcf,'-property','FontSize'),'FontSize',20,'fontweight','bold')

toc

function [yprime] = NOCmodel2(t,y,T)
Rm = y(1);
N  = y(2);

%A. Known variables
Xp         = 1; 
Xo         = 0.06; 
R          = 8.314 ;                             %Gas constant
kB         = 1.38064852 * 1e-23 ;                %R/avogadro no.
Na         = 6.02214 * 1e23;                     %avogadro no.
IncXe      = 0.01;                               %Equilibrium mol frac.
sigma      = 0.13;                               %interfacial free energy J/m2
Di         = 5e-20;                              %aDiffusivity
Vinc       = 1.6*1e-29*Na;                       %m3/mol
VFe        = Vinc;                               %Vinc;


%B. Calculated properties

alpha      = VFe/Vinc;
Rzero      = 2*sigma*Vinc/R/T;
IncXieq(1) = IncXe ;
IncXi(1)   = IncXieq(1)*exp(2*sigma*Vinc/(R*Rm*T));

if Rm <0
   error('negative Rm') 
elseif N < 0
   error('negative N') 
end

% Solute balance
X_O = (Xo - alpha*(4/3)*pi*Rm^3*N*Xp)/(1 - (alpha*(4/3)*pi*Rm^3*N*Xp));


%C. number density
[~,~,Rstar,~,Rcrit,~,C1] = nucleate2(X_O, T, t) ;

f_coars = 0 ;
if Rm < 1.01*Rstar && Rm > 0.99*Rstar
    f_coars = 1-(1000*((Rm/Rstar) - 1)^2) ;
end

D1 = (4/27) * ((IncXieq(1))/(alpha*Xp - IncXieq(1))) * (Rzero*Di/Rm^3) * ((((Rzero*IncXieq(1))/(Rm*(Xp - IncXieq(1)))) * ((3/(4*pi*Rm^3))) - N) - 3*N);

dNdt = abs(min(0,(C1+D1)/abs(C1+D1)))*f_coars*D1 + max(0,(C1+D1)/abs(C1+D1))*C1;


%D. mean radius
A1 = (Di/Rm)*((X_O - IncXi(1))/(alpha*Xp - IncXi(1))); A2 = (1/N)*dNdt*(Rcrit - Rm);

A = A1 + A2 ;
B = (4/27) * ((IncXieq(1))/(alpha*Xp - IncXieq(1))) * ((Rzero)*Di/Rm^2) ;

%fprintf('t = %1.1e, A1 =  %2.2e, A2 =  %2.2e, B =  %2.2e, C =  %2.2e, D =  %2.2e, fc = %2.2f, Rm = %2.2e, dNdt = %2.2e, xo = %2.2e, N = %2.2e, Rcrit = %2.2e\n',t,A1,A2,B,C1,D1,f_coars,Rm, dNdt, X_O,N,Rcrit )

yprime = [f_coars*B + (1-f_coars)*A ; abs(min(0,(C1+D1)/abs(C1+D1)))*f_coars*D1 + max(0,(C1+D1)/abs(C1+D1))* C1];


end

function [SS,delGv,rc,Z,rstar,delGstar,I_nucl] = nucleate2(Xo, T, t)
R         = 8.314 ;                             %Gas constant
kB        = 1.38064852 * 1e-23 ;                %R/avogadro no.
Na        = 6.02214 * 1e23;                     %avogadro no.
IncXe     = 0.01;                               %Equilibrium mol frac.
sigma     = 0.13;                               %interfacial free energy J/m2
Di        = 5e-20;                              %aDiffusivity
Vinc      = 1.6*1e-29*Na;
VFe       = Vinc;                               %Vinc;
a         = 4.04e-10 ;                          %lattice parameter FCC iron 


Vm1       =  Vinc;                              %m3/mol  
Vm        =  Vm1/Na;                            %m3
SS        =  Xo/IncXe;                          %reaction quotient of bulk steel
delGv     = -R*T*log(SS)/Vm1;                   %chemical driving force (J/m3)                                                                   
rc        = -2*sigma/delGv ;                    %critical radius (m)
Z         =  Vm/(2*pi*rc^2)*sqrt(sigma/kB/T);   %Zeldovich factor
rstar     =  rc + 0.5*sqrt(kB*T/(pi*sigma)) ;   %(1/(Z*sqrt(pi))); %m
delGstar  =  16*pi*(sigma)^3/(3*delGv^2) ;      % J



betastar  = 4*pi*rc^2*Di*Xo/a^4;
tow       = 0.5/(betastar*Z); 

PreEx     = betastar * Z * Na/VFe;
I_nucl    = PreEx * exp(-delGstar/kB/T) * exp(-tow/t);



end

}

$\endgroup$
14
  • 1
    $\begingroup$ If you could add the equations to your question, you might get an answer that's more informative then "It all depends on the details of your problem." $\endgroup$ – Brian Borchers Feb 23 at 20:12
  • 1
    $\begingroup$ In theory it should not, but in practice the solvers make heuristical decisions generalizing the behavior of test cases. So if the non-dimensionalizing moves your system towards these test cases, which it should if done sensibly, then the behavior of the solver can improve, which might look like stiffness was reduced. $\endgroup$ – Lutz Lehmann Feb 23 at 21:23
  • $\begingroup$ Hi, I have added the eqns for refernce. @LutzLehmann, thanks for the clarification. $\endgroup$ – Ronnie1993 Feb 23 at 22:32
  • 2
    $\begingroup$ It looks like $N'$ has a jump discontinuity in its ODE. How do you treat this? Branching in the ODE is sub-optimal, it is better to have two smooth ODE systems and use an event mechanism to switch. $\endgroup$ – Lutz Lehmann Feb 23 at 22:45
  • 1
    $\begingroup$ Yes, that produces a jump. The solver "sees" the jump as a region of stiff behavior and thus reduces the step size. If you enter a sliding mode, then this persists for some time. If you manage the switch externally, this disappears. If you are careful, you can get a similar improvement by mollifying the jump into some continuous function with a steep slope. See stackoverflow.com/questions/60309851, stackoverflow.com/questions/54767096 for examples $\endgroup$ – Lutz Lehmann Feb 24 at 7:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.