Get wrong answers with 2D Newton iteration method

I need some help with the following problem involving the Newton method.

Basically I tried to use a model called LGK model to predict the dendrite growth velocity and dendrite tip radius during the solidification of an alloy. First I wanted to use the same parameters that the model creators used, to see if I can get the same results as theirs. A 2D Newton method is used to iterate until the result is reached.

However, despite I used exactly the same parameters and equations (I have checked them multiple times), I could not get the same result as the model creators, instead I got NAN after several cycles.

Here is my code:

clear
clc

%Parameters:
%specific heat capacity (J/kg*K)
Cp = 1937.5;
%latent heat (J/kg)
deltaH = 46.26*10^3;
%thermal diffusivity (m^2/s)
a = 1.14*10^-7;
%Gibbs-Thomson coefficient (K*m)
gamma = 6.62*10^-8;
%interdiffusion coefficient (m^2/s)
D = 1.27*10^-9;
%liquid slope (K/mol.%)
m = -2.16;
%partition coefficient
k0 = 0.103;
%stability constant
sigma = 1/(4*(pi^2));

%composition and temperature
C0 = 0;
deltaT = 0.5;

%initial values to start iterating
V = 0.0001;
R = 0.00001;

syms x1 x2
Pt = (x1*x2)/(2*a); %thermal Peclet number
Pc = (x1*x2)/(2*D); %solutal Peclet number
IvPc = Pc*exp(Pc)*expint(Pc); %Ivantsov function for solutal field
IvPt = Pt*exp(Pt)*expint(Pt); %Ivantsov function for thermal field
f1 = (deltaH/Cp)*IvPt - m*C0*((1-k0)*IvPc)/(1-(1-k0)*IvPc) + 2*gamma/x2 - deltaT;
f2 = (gamma/sigma)/((Pt*deltaH/Cp)-Pc*m*C0*(1-k0)/(1-(1-k0)*IvPc)) - x2;

df1x1 = diff(f1,x1);
df1x2 = diff(f1,x2);
df2x1 = diff(f2,x1);
df2x2 = diff(f2,x2);

%iteration
f3 = [x1;x2] - (([df1x1 df1x2; df2x1 df2x2])^-1)*[f1;f2];
for n = 1:20
x1 = V; x2 = R;
ans_in_sym = subs(f3);
ans_in_double = double(ans_in_sym);
V = ans_in_double(1);
R = ans_in_double(2);
end

V_Difference = V-x1;
R_Difference = R-x2;
fprintf('V = %e\n R = %e\n V-V0 = %e\n R-R0 = %e',V,R,V_Difference,R_Difference);

For these specific C0 and deltaT, V and R should be around 4*10^-5 and 26*10^-6 respectively, but with my code I got NAN after about 45 cycles. First I thought the structure of iteration might be wrong, so I tried a much easier equation:

f1 = 3*x1^2+2*x2^2-25;
f2 = 2*x1^2-x2-15;

and I successfully got the answer with the same (except the equation part) code. So I think the iteration structure is fine, and the problem should be in the equation part.

But I have checked the equation multiple times... Is it possible that diff() gives wrong partial differential expression for a complicated symbolic expression? I think that might be the problem because I got pretty long symbolic expressions in my code.

Thank you very much!

Answer 1

[Benjamin Thompson]

There are a lot of really large numbers in these functions, leading to very steep gradient calculations. It may just not be numerically stable. Can you pick units so that parameters do not need such large exponents?