Solving a transcendental equation

Hello,

I am trying to solve a transcendental equation for solidification process. You can see the equation below.

All parameters are known except "lambda" and I need to solve this equation to calculate "lambda". The problem is that I cannot be sure how to solve this equation by using MATLAB. I have tried the code given below:

%% Known Parameters

% Pyhsical Parameters
L = 0.02; %[m]
H = 0.5; %[m]

%Thermophysical Parameters
T_i = 10; %[C]
T_0 = -5; %[C]
T_m = 0; %[C]

c_solid = 2050; %[J.kg-1.K-1]
c_liquid = 4181.3; %[J.kg-1.K-1]

rho_solid = 1000; %[kg/m^3]
rho_liquid = 999.7026; %[kg/m^3]

k_solid = 0.5610; %[W/m.K]
k_liquid = 0.5674; %[W/m.K]

%% Calculation of Constants

alpha_solid = k_solid / (rho_solid * c_solid);
alpha_liquid = k_liquid / (rho_liquid * c_liquid);

%% Solver
syms lambda

eqn = exp(-lambda^2) / erf(lambda) + (k_liquid / k_solid) * sqrt((alpha_solid / alpha_liquid)) * ((T_m - T_i) / (T_m - T_0)) * exp(-lambda^2 * (alpha_solid / alpha_liquid)) / erfc(lambda * sqrt((alpha_solid / alpha_liquid))) == lambda * L * sqrt(pi) / (c_solid * (T_m - T_0));

solve(eqn, lambda)

This code returns a negative answer which is wrong. I beleive I am using a wrong function to solve the equation. If you could help me it would be very appreciated.

Regards

Answer 1

[David Goodmanson]

Hi Kubilay,

Finding an algebraic solution for this equation in terms of lambda is not going to happen. But a numerical solution works fine:

% look for a solution that you can see
lambda = 0:.01:1;
plot(lambda,fun(lambda))
lambda0 = fzero(@fun,[.1 .3])
function y = fun(lambda)
% Pyhsical Parameters
L = 0.02; %[m]
H = 0.5; %[m]
%Thermophysical Parameters
T_i = 10; %[C]
T_0 = -5; %[C]
T_m = 0; %[C]
c_solid = 2050; %[J.kg-1.K-1]
c_liquid = 4181.3; %[J.kg-1.K-1]
rho_solid = 1000; %[kg/m^3]
rho_liquid = 999.7026; %[kg/m^3]
k_solid = 0.5610; %[W/m.K]
k_liquid = 0.5674; %[W/m.K]
% Calculation of Constants
alpha_solid = k_solid / (rho_solid * c_solid);
alpha_liquid = k_liquid / (rho_liquid * c_liquid);

y1 = exp(-lambda.^2) ./ erf(lambda) + (k_liquid / k_solid) ...
     * sqrt((alpha_solid / alpha_liquid)) * ((T_m - T_i) / (T_m - T_0)) ... 
     * exp(-lambda.^2 * (alpha_solid / alpha_liquid))...
     ./erfc(lambda* sqrt((alpha_solid / alpha_liquid)));
y2 = lambda * L * sqrt(pi) / (c_solid * (T_m - T_0));
y = y1-y2;
end
   lambda0 =
     0.2177

[Paul]

Hi David,

A slight modification that allows the user to define eqn symbolically, perhaps for additional analysis, but still solve using fzero

% Pyhsical Parameters
L = 0.02; %[m]
H = 0.5; %[m]
%Thermophysical Parameters
T_i = 10; %[C]
T_0 = -5; %[C]
T_m = 0; %[C]
c_solid = 2050; %[J.kg-1.K-1]
c_liquid = 4181.3; %[J.kg-1.K-1]
rho_solid = 1000; %[kg/m^3]
rho_liquid = 999.7026; %[kg/m^3]
k_solid = 0.5610; %[W/m.K]
k_liquid = 0.5674; %[W/m.K]
%% Calculation of Constants
alpha_solid = k_solid / (rho_solid * c_solid);
alpha_liquid = k_liquid / (rho_liquid * c_liquid);
%% Solver
syms lambda
eqn = exp(-lambda^2) / erf(lambda) + (k_liquid / k_solid) * sqrt((alpha_solid / alpha_liquid)) * ((T_m - T_i) / (T_m - T_0)) * exp(-lambda^2 * (alpha_solid / alpha_liquid)) / erfc(lambda * sqrt((alpha_solid / alpha_liquid))) == lambda * L * sqrt(pi) / (c_solid * (T_m - T_0));
eqn2solve = matlabFunction(lhs(eqn)-rhs(eqn));
fplot(eqn2solve,[0 3]);grid

fzero(eqn2solve,[.1 .3])

ans = 0.2177

[Kubilay AKPINAR]

Dear David and Dear Paul, thank you for your explanation. Everything is more clear now. Appreciated.

Answer 2

[Paul]

Here's another way to get a numeric solution, staying in the symbolic world

% Pyhsical Parameters
L = 0.02; %[m]
H = 0.5; %[m]
%Thermophysical Parameters
T_i = 10; %[C]
T_0 = -5; %[C]
T_m = 0; %[C]
c_solid = 2050; %[J.kg-1.K-1]
c_liquid = 4181.3; %[J.kg-1.K-1]
rho_solid = 1000; %[kg/m^3]
rho_liquid = 999.7026; %[kg/m^3]
k_solid = 0.5610; %[W/m.K]
k_liquid = 0.5674; %[W/m.K]
%% Calculation of Constants
alpha_solid = k_solid / (rho_solid * c_solid);
alpha_liquid = k_liquid / (rho_liquid * c_liquid);
%% Solver
syms lambda
eqn = exp(-lambda^2) / erf(lambda) + (k_liquid / k_solid) * sqrt((alpha_solid / alpha_liquid)) * ((T_m - T_i) / (T_m - T_0)) * exp(-lambda^2 * (alpha_solid / alpha_liquid)) / erfc(lambda * sqrt((alpha_solid / alpha_liquid))) == lambda * L * sqrt(pi) / (c_solid * (T_m - T_0));

Plot the equation, get an idea where the roots might be:

fplot(lhs(eqn)-rhs(eqn),[0 10])

Looks like the range of [0 1] should work

vpasolve(eqn,lambda,[0 1])

ans = 

0.21772684787123896059255753602836