Solving a transcendental equation

164 次查看(过去 30 天)
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

采纳的回答

David Goodmanson
David Goodmanson 2021-11-3
编辑:David Goodmanson 2021-11-3
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 个评论
Paul
Paul 2021-11-4
编辑:Paul 2021-11-4
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 2021-11-4
Dear David and Dear Paul, thank you for your explanation. Everything is more clear now. Appreciated.

请先登录,再进行评论。

更多回答(1 个)

Paul
Paul 2021-11-4
编辑:Paul 2021-11-4
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

类别

Help CenterFile Exchange 中查找有关 Video games 的更多信息

标签

产品


版本

R2020a

Community Treasure Hunt

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

Start Hunting!

Translated by