# Errors using function and fsolve

27 views (last 30 days)
Brady Wayne Robinson on 24 Jan 2020 at 4:07
Commented: Walter Roberson on 24 Jan 2020 at 5:43
Below is the code I have written and then I am getting the error messages. I am assuming my problem lies within the function or the fsolve. I have tried and tried and tried and cannot figure out what is wrong. Thank You for help!
I have tried so hard and would really be grateful for help! Thank You so much!
% *Set universal constants:*
sigma = 5.67e-8; % Stefan-Boltzmann constant (in W/m^2/K^4)
%
% *Set surrounding properties*
surr.T_flow = 295; % surrounding flow temperature (in K)
surr.T_wall = 295; % surrounding wall temperature (in K)
%
% *Set propeties for the computer chip*
chip.dy = 0.005; % chip vertical thickness (in m)
chip.A_surf = 0.002; % chip surface area in contact with heat sink (in m^2)
chip.k = 20; % chip thermal conductivity (in W/m/K)
%
% *Set propeties for the heat sink*
sink.eps = 0.7; % heat sink
sink.A_surf = 0.02; % heat sink surface area (in m^2)
%%
% *Make vectors of convection coefficients and chip power*
surr.h =[6,11,16,21,26] % heat transfer coefficient (in W/m^2*K)
n_h = length(surr.h);
P_elec = [0:1:30] % chip power dissipation (in W)
n_P = length(P_elec);
sink.T = zeros(n_P,n_h); % allocate memory for sink temperatures (in K)
for i_h=1:n_h
for i_P=1:n_P
%%
% *Call the residual function to solve for sink temperature using f_zero*
T_guess=300
sink.T(i_h) = fsolve(@(T) Resfun(T,surr.h(i_h),sink,surr,sigma),T_guess);
sink.T(i_P) = fsolve(@(T) Resfun(T,P_elec(i_P),sink,surr,sigma),T_guess);
end
end
% *Residual function for finding sink temperature*
function res = Resfun(T,P_elec,sink,surr,sigma)%inputs
% set up energy balance residual
res= P_elec - surr.h.*sink.A_surf.*(T-surr.T_flow)-sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4); % energy balance eq
end

#### 1 Comment

Brady Wayne Robinson on 24 Jan 2020 at 4:26
There it is! Thank You

Walter Roberson on 24 Jan 2020 at 4:15
Edited: Walter Roberson on 24 Jan 2020 at 4:16
The function you are trying to find the zero of, must return a scalar. Your formula involves surr.h which is a vector, so your function is returning a vector.
Note that fzero which you used has this restriction, but fsolve does not have this restriction.

Brady Wayne Robinson on 24 Jan 2020 at 4:26
I tried using Fsolve and still got a boatload of errors
Walter Roberson on 24 Jan 2020 at 5:14
If you have the symbolic toolbox then
syms T
F = Resfun(T,P_elec(i_P),sink,surr,sigma)
S = solve(F(1),T);
vpa(subs(F(2:end), T, S))
You will see a matrix of non-zero values, too large in value to be just round off error. This tells you that the solution to the first of your equations cannot satisfy the other four equations. Therefore fsolve() is not going to be able to find a T value that solves all 5 equations simultaneously.
Your equations are hopeless in the form written.
Walter Roberson on 24 Jan 2020 at 5:27
Might I suggest to you that you do not want to try to find a T that satisfies all
surr.h =[6,11,16,21,26] % heat transfer coefficient (in W/m^2*K)
simultaneously, and that instead you want to loop through solving your system for each one of those individually.
You seem to be giving that a try with
n_h = length(surr.h);
for i_h=1:n_h
...
sink.T(i_h) = fsolve(@(T) Resfun(T,surr.h(i_h),sink,surr,sigma),T_guess);
so you are passing in the individual surr.h value as the second parameter. But inside Resfun you have
res= P_elec - surr.h.*sink.A_surf.*(T-surr.T_flow)-sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4); % energy balance eq
which uses all of surr.h
By the way,
sink.T(i_h) = fsolve(@(T) Resfun(T,surr.h(i_h),sink,surr,sigma),T_guess);
does not involve variables from the inner loop for i_P=1:n_P so there is no obvious point in doing inside the inner loop.
And
sink.T(i_P) = fsolve(@(T) Resfun(T,P_elec(i_P),sink,surr,sigma),T_guess);
is going to overwrite the same location at some point.
Perhaps you should be saving the results to a 2D array indexed at i_h and i_P . But you do need to resolve the question of why you are calling the function with two different inputs and writing to the same output variable.

Ridwan Alam on 24 Jan 2020 at 5:27
% *Set universal constants:*
global surr sink sigma
sigma = 5.67e-8; % Stefan-Boltzmann constant (in W/m^2/K^4)
%
% *Set surrounding properties*
surr.T_flow = 295; % surrounding flow temperature (in K)
surr.T_wall = 295; % surrounding wall temperature (in K)
%
% *Set propeties for the computer chip*
chip.dy = 0.005; % chip vertical thickness (in m)
chip.A_surf = 0.002; % chip surface area in contact with heat sink (in m^2)
chip.k = 20; % chip thermal conductivity (in W/m/K)
%
% *Set propeties for the heat sink*
sink.eps = 0.7; % heat sink
sink.A_surf = 0.02; % heat sink surface area (in m^2)
% *Make vectors of convection coefficients and chip power*
surr.h =[6,11,16,21,26]; % heat transfer coefficient (in W/m^2*K)
n_h = length(surr.h);
P_elec = [0:1:30]; % chip power dissipation (in W)
n_P = length(P_elec);
sink.T = zeros(n_P,n_h); % allocate memory for sink temperatures (in K)
for i_h=1:n_h
for i_P=1:n_P
% *Call the residual function to solve for sink temperature using f_zero*
T_guess=300;
temp_var = fsolve(@Resfun,[T_guess,surr.h(i_h)])
sink.T(i_P,i_h) = temp_var(1);
temp_var = fsolve(@Resfun,[T_guess,P_elec(i_P)])
sink.T(i_P,i_h) = temp_var(1);
end
end
% *Residual function for finding sink temperature*
function res = Resfun(alpha)%inputs
global surr sink sigma
T = alpha(1);
P_elec_i = alpha(2);
%sink = alpha{3};
%surr = alpha{4};
%sigma = alpha{5};
% set up energy balance residual
res= P_elec_i - ...
surr.h.*sink.A_surf.*(T-surr.T_flow) - ...
sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4);
% energy balance eq
end
Hope this helps.

Walter Roberson on 24 Jan 2020 at 5:40
global almost always causes more problems than it solves.
Walter Roberson on 24 Jan 2020 at 5:43
function res = Resfun(alpha)%inputs
global surr sink sigma
T = alpha(1);
P_elec_i = alpha(2);
Is the goal now to have Resfun vary T and P_elec_i to try to find a zero of the equation?
If so then calculate
surr.h.*sink.A_surf.*(T-surr.T_flow) - ...
sink.A_surf.*sigma.*sink.eps.*(T^4-surr.T_wall^4)
and set P_elec_i to that value, and you get a perfect 0 output.
Well, that is, you would, if you had fixed the problem that you are calculating for all surr.h values simultaneously, and I can guarantee you that you will not find a single [T, P_elect_i] pair that will solve all 5 equations simultaneously.