Solving a transcendental equation

28 views (last 30 days)
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

Accepted Answer

David Goodmanson
David Goodmanson on 3 Nov 2021
Edited: David Goodmanson on 3 Nov 2021
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
  2 Comments
Kubilay AKPINAR
Kubilay AKPINAR on 4 Nov 2021
Dear David and Dear Paul, thank you for your explanation. Everything is more clear now. Appreciated.

Sign in to comment.

More Answers (1)

Paul
Paul on 4 Nov 2021
Edited: Paul on 4 Nov 2021
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

Products


Release

R2020a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!