Solve a function and plot its contour plot. Not getting the desired contour plot?

I have defined the temperature field as Z..and want to plot the temperature contour. However, I am unable to get the desired contour plot. Can someone please help me with this? I have also trield fcontour by defining X,Y as variables..but with no results.

P = 50;
v = 0.1;
k = 113;
Tm = 843;
T0 = 300;
a = 4.63 * 10^(-5);
eps = 0.9;
sig = 5.67 * 10^-8;
A = 10^-5;
kp = 0.21;
x = linspace(-3, 3);
y = linspace(-3, 0);
% Remove NaN values by replacing them with a default value (e.g., 0)
x(x == 0) = 0;
y(y == 0) = 0;
[X, Y] = meshgrid(x, y);
% Ensure that r is not zero to avoid division by zero issues
r = sqrt((X.*(10^-3)).^2 + (Y.*(10^-3)).^2);
r(r == 0) = 10^-6; % Replace zeros with a small value (eps) to avoid division by zero
Z = (1./(4*k*pi.*r.*(Tm-T0))) .* (P * exp((-v.*(r+X.*10^-3))./(2*a)) - A*(h.*(Z-T0)+eps*sig*(Z.^4-T0^4)+(kp.*(Z-T0)./r)));
figure
contourf(X, Y, Z)
colorbar;

[Dyuman Joshi]

What is the desired output?

h has not been defined for the above code.

Additionally, eps is a built-in function in MATLAB. It's better not to name variables (or scripts for that matter) using functions. I suggest you rename that variable.

% Remove NaN values by replacing them with a default value (e.g., 0)
x(x == 0) = 0;
y(y == 0) = 0;

This is redundant. Why have you used this? And that code does not do what is written the comment above it. So, it is confusing to me.

[Voss]

h = 1;
Z = 1000;

P = 50;
v = 0.1;
k = 113;
Tm = 843;
T0 = 300;
a = 4.63 * 10^(-5);
eps = 0.9;
sig = 5.67 * 10^-8;
A = 10^-5;
kp = 0.21;
x = linspace(-3, 3);
y = linspace(-3, 0);
% Remove NaN values by replacing them with a default value (e.g., 0)
x(x == 0) = 0;
y(y == 0) = 0;
[X, Y] = meshgrid(x, y);
% Ensure that r is not zero to avoid division by zero issues
r = sqrt((X.*(10^-3)).^2 + (Y.*(10^-3)).^2);
r(r == 0) = 10^-6; % Replace zeros with a small value (eps) to avoid division by zero
Z = (1./(4*k*pi.*r.*(Tm-T0))) .* (P * exp((-v.*(r+X.*10^-3))./(2*a)) - A*(h.*(Z-T0)+eps*sig*(Z.^4-T0^4)+(kp.*(Z-T0)./r)));
figure
contourf(X, Y, Z)
colorbar;

[Walter Roberson]

% Remove NaN values by replacing them with a default value (e.g., 0)
x(x == 0) = 0;
y(y == 0) = 0;

If you want to remove NaN, here are some ways:

%way 1
x(x ~= x) = 0;
y(y ~= y) = 0;

%way 2
x(isnan(x)) = 0;
y(isnan(y)) = 0;

%way 3
x = fillmissing(x, 'Constant', 0);
y = fillmissing(y, 'Constant', 0);

%way 4, more robust
mask = isnan(x) | isnan(y);
x(mask) = 0;
y(mask) = 0;

Note, however, that x = linspace(-3, 3); will not generate any NaN values. You should probably be post-processing Z instead.

[AD]

@Voss Thanks for helping! But why have you initialized the variable Z initially?

[Walter Roberson]

Your code has

Z = (1./(4*k*pi.*r.*(Tm-T0))) .* (P * exp((-v.*(r+X.*10^-3))./(2*a)) - A*(h.*(Z-T0)+eps*sig*(Z.^4-T0^4)+(kp.*(Z-T0)./r)));

You have the problem that you have not defined Z but you use Z on the right hand side.

Are you trying to say that you want to calculate Z such that there is equality between Z and the right hand side calculation involving Z?

[AD]

Yes..I want to solve for Z for each (X,Y) and plot its contour plot.

[Walter Roberson]

Because of the Z.^4 on the right hand size, you are defining a quartic -- a polynomial in degree 4. There are 4 solutions for each point. An even number of those solutions will be real-valued.

[AD]

So, is it possible to plot only one of the real valued solution?

Answer 1

[Walter Roberson]

%h was not defined in original code -- make sure you assign a meaningful
%value!
h = 1;


syms X Y Z real
Q = @(v) sym(v);

P = Q(50);
v = Q(0.1);
k = Q(113);
Tm = Q(843);
T0 = Q(300);
a = Q(463) * Q(10)^(-7);
eps = Q(0.9);
sig = Q(567) * Q(10)^-10;
A = Q(10)^-5;
kp = Q(0.21);
Pi = Q(pi);

R = sqrt((X.*(Q(10)^-3)).^2 + (Y.*(Q(10)^-3)).^2);
r = piecewise(R == 0, 1e-6, R);

eqn = Z == (1./(4*k*Pi.*r.*(Tm-T0))) .* (P * exp((-v.*(r+X.*10^-3))./(2*a)) - A*(h.*(Z-T0)+eps*sig*(Z.^4-T0^4)+(kp.*(Z-T0)./r)));
zsol = solve(eqn, Z, 'returnconditions', true)

zsol = struct with fields:
             Z: [6×1 sym]
    parameters: [1×0 sym]
    conditions: [6×1 sym]

x = linspace(-3, 3);
y = linspace(-3, 0);
[xG, yG] = meshgrid(x, y);

%warning: zsolfun returns a matrix and must be invoked on scalars!
zsolfun = matlabFunction(reshape(zsol.Z, 1, 1,[]), 'File', 'zsol.m', 'Vars', [X, Y], 'optimize', false);
zcondfun = matlabFunction(reshape(zsol.conditions, 1, 1, []), 'File', 'zcond.m', 'Vars', [X, Y], 'optimize', false);
[xG, yG] = meshgrid(x, y);
Zcell = arrayfun(zsolfun, xG, yG, 'uniform', 0);
Zmat = cell2mat(Zcell);
Zcondcell = arrayfun(zcondfun, xG, yG, 'uniform', 0);
Zcond = cell2mat(Zcondcell);

for L = 1 : size(Zcond,3)
    mask = ~Zcond(:,:,L);
    layer = Zmat(:,:,L);
    layer(mask) = NaN;
    if nnz(~isnan(layer)) == 0; continue; end
    figure;
    subplot(2,1,1)
    contour(xG, yG, layer, 7);
    colorbar();
    title("root #" + L);
    subplot(2,1,2)
    scatter(xG(:), yG(:), [], layer(:));
    colorbar();
end