无量纲 ODE 会影响方程的刚度吗？能否提高MATLAB中ode45、ode113等数值方法的稳定性？
我正在尝试解决 2 个 eqns。这可能需要非常小的步长才能获得像 1e-11 这样的稳定性。我想提高数值方法的速度。非维度化在这方面有帮助吗？
ODE（IV：Rbar，N）显示在链接1中。输出应接近2所示
我的测试问题的工作代码：
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
}