Solving a transcendental equation

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

Categories

Find more on Video games in Help Center and File Exchange

Products


Release

R2020a

Community Treasure Hunt

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

Start Hunting!