How can I write and solve this equation? as a function??

[Benjamin Thompson]

Solve the equation for what? Need more information, and the screenshot is a little difficult to read.

Answer 1

[Star Strider]

Possibly, after first providing values for the constants:

T(z,z0,Tz0,Tinf) = Tinf - (Tinf - Tz0).*exp(z-z0);              % Function
zsol = fzero(@(z)T(z,z0,Tz0,Tinf), rand)                        % Correct Call To It (Here Using 'fzero', Can Be Any Function)

See the documentation on Anonymous Functions for details.

.

[Cesar Cardenas]

Thanks much for your help. I need to calculate the temperature (T) at 200 (z), z0 is 120, Tinf is 651 and Tz0 is 350. So, in fzero, woud I have to type f(200)?? not sure about that??.

Additionally, how I would have to type this equation to compute s especically dT/dz?? as a derivative?

thanks much.

[Walter Roberson]

T(z,z0,Tz0,Tinf) = Tinf - (Tinf - Tz0).*exp(z-z0);              % Function

That syntax would require Symbolic Toolbox and would define a Symbolic Function

zsol = fzero(@(z)T(z,z0,Tz0,Tinf), rand)                        % Correct Call To It (Here Using 'fzero', Can Be Any Function)

Calls to Symbolic Functions return symbolic values, but fzero() needs the value to be floating point. You would need

zsol = fzero(@(z)double(T(z,z0,Tz0,Tinf)), rand)

[Walter Roberson]

format long g
Tinf = 651; Tz0 = 350; z0 = 20;
Tz = @(z) Tinf - (Tinf - Tz0).*exp(z-z0);
Tz(200)

ans = 
     -4.48304644435333e+80

[Star Strider]

As always, my pleasure!

The function I originally wrote was for the Symbolic Math Toolbox, although the fzero call was numeric.

The numeric version (written as such) produces some interesting results —

Tfcn = @(z,z0,Tz0,Tinf) Tinf - (Tinf - Tz0).*exp(z-z0);            % Function

T = 200;
z0 = 120;
Tinf = 651;
Tz0 = 350;

[zsol,fv] = fzero(@(z)(T - Tfcn(z,z0,Tz0,Tinf)), 100)               % Correct Call To It (Here Using 'fzero', Can Be Any Function)

zsol = 120.4044

fv = 6.2528e-13

[zsol,fv] = fsolve(@(z)(T - Tfcn(z,z0,Tz0,Tinf)), 100)              % Correct Call To It (Here Using 'fzero', Can Be Any Function)

Equation solved.

fsolve completed because the vector of function values is near zero
as measured by the value of the function tolerance, and
the problem appears regular as measured by the gradient.

zsol = 120.4044

fv = 6.2528e-13

zv = linspace(zsol - 1.5, zsol + 1.5, 500);
zi = interp1(Tfcn(zv,z0,Tz0,Tinf), zv, 0)

zi = 120.7714

figure
plot(zv, Tfcn(zv,z0,Tz0,Tinf))
hold on
plot(zsol, fv, 'cp', 'MarkerSize',10)
plot(zi, fv, 'mp', 'MarkerSize',10)
hold off
grid
xlabel('z')
ylabel('T(z)')
legend('Function Evaluation','Function Solution','Interpolation Solution', 'Location','best')

I leave you to explore the solution you prefer.

It might also be interesting to see what the Symbolic Math Toolbox solution is.

(Computer crashed — again — a few minutes ago, explaining much of the delay. I have no idea the reason both Windows 10 and Windows 11 do this so frequently.)

.

[Walter Roberson]

dTdz = @(z) appropriate calculation of derivative
s = dTdz(z0) / (Tinf - Tz0)

[Star Strider]

@Walter Roberson — Thank you!

I forgot about the derivative — it has just been one of those days.

The approach I usually use is a sort of central difference approximation —

T = 200;
z0 = 120;
Tinf = 651;
Tz0 = 350;

Tfcn = @(z,z0,Tz0,Tinf) Tinf - (Tinf - Tz0).*exp(z-z0);            % Function
h = 1E-8;
dTdz = @(z) (Tfcn(z+h,z0,Tz0,Tinf) - Tfcn(z-h,z0,Tz0,Tinf)) / (2*h);
z = 120.7714;
Derivative = dTdz(z)

Derivative = -651.0000

zv = linspace(z - 1.5, z + 1.5, 500);
Check = interp1(zv, gradient(Tfcn(zv,z0,Tz0,Tinf))./gradient(zv), z)

Check = -651.0073

.

[Cesar Cardenas]

Thanks much for your help, a variable is missing (-s) in this equation ( as shown in the image): so it would be

Tfcn = @(z,z0,Tz0,Tinf) Tinf - (Tinf - Tz0).*exp-s(z-z0);

So, I'm including -s and gave a value, but I get this error

Error in upperatmospherehw2 (line 8)

[zsol,fv] = fzero(@(z)(T - Tfcn(z,z0,s,Tz0,Tinf)), 100)

Tfcn = @(z,z0,s,Tz0,Tinf) Tinf - (Tinf - Tz0).*exp-s(z-z0);            % Function
T = 200;
z0 = 120;
Tinf = 651;
Tz0 = 350;
s = 0.5;

[zsol,fv] = fzero(@(z)(T - Tfcn(z,z0,s,Tz0,Tinf)), 100)    

[Cesar Cardenas]

For this equation I have:

Tz0 = 350;
z0 = 120;
Tinf = 651;

I just need to calculate s, normally dTdz is a small value.

[Cesar Cardenas]

Thanks much for your help. This equation is missing the variable s As shown in the image). I added and assigned a value. However, I'm getting this error not sure why??

Error in upperatmospherehw2 (line 29)

Tz(200)

format long g
Tinf = 651; Tz0 = 350; z0 = 20; s = 0.021
Tz = @(z) Tinf - (Tinf - Tz0).*exp-s*(z-z0);
Tz(200)

[Star Strider]

I decided to switch to the Symbolic Math Toolbox because of the derivative.

The image was so small that I missed the ‘s’.

Of interest, MATLAB gets a different result with the integration —

syms s T(z) T_inf T_z0 z z_0 
D1T = diff(T);
Eqn = D1T == s*(T_inf - T_z0)

Eqn(z) = 

Tf(z) = simplify(dsolve(Eqn), 500)

Tf(z) = 

There is something about that equation that is not obvious here.

That aside, solving for ‘s’:

clear vars
syms s Tz T_inf T_z0 z z_0
Eqn = Tz == T_inf - (T_inf - T_z0) * exp(-s*(z-z_0))

Eqn = 

D1Tz = diff(rhs(Eqn),z)

D1Tz = 

Eqn = isolate(Eqn, s)

Eqn = 

Tz = 200;
z_0 = 120;
T_inf = 651;
T_z0 = 350;
z = 120.7714;
Eqn = subs(Eqn)

Eqn = 

Eqn = vpa(Eqn)

Eqn = 

s = rhs(Eqn)

s = 

D1Tz = subs(D1Tz)

D1Tz = 

D1Tz = vpa(D1Tz)

D1Tz = 

This does not look small to me, however perhaps in the context of this system, it is.

.

[Cesar Cardenas]

Thanks so much, however, I think that the value for z that you got earlier (z = 120.7714) would not be right because you missed ¨"s". I´ve been trying to fix it but still get an error, not sure why??

Tz = T_inf - (T_inf - T_z0) * exp(-s*(z-z_0))           
T = 200;
z0 = 120;
Tinf = 651;
Tz0 = 350;
s = 0.021;

[zsol,fv] = fzero(@(z)(T - Tfcn(z,z0,Tz0,Tinf,s)), 100)  

[Star Strider]

The ‘Tfcn’ code needs to change (and ‘Tfcn’ was not in the code for your most recent Comment, so that could have been the error).

This appear to work, and gives the same result for several different initial estimates —

Tfcn = @(z,z0,s,Tz0,Tinf) Tinf - (Tinf - Tz0).*exp(-s.*(z-z0));            % Function
T = 200;
z0 = 120;
Tinf = 651;
Tz0 = 350;
% s = 0.5;
s = 0.021;

[zsol,fv] = fzero(@(z)(T - Tfcn(z,z0,s,Tz0,Tinf)), 150)

zsol = 100.7449

fv = 1.1369e-13

Tcalc = Tfcn(zsol,z0,s,Tz0,Tinf)

Tcalc = 200.0000

I assume the curve is similar to the previous plot, so I will not re-plot it here.

.

[Cesar Cardenas]

great thanks so much, it works well... how can I visualize or see the values you are getting for s? I'm trying to get the value of s for z0 = 110...not sure why it is negative since s would be basically a distance??

[Star Strider]

As always, my pleasure!

I am not certain what it would mean to visualilse them, however one way would be to choose one or two specific variables as independent variables, create vectors for one or both them (perhaps ‘z’ and ‘T’, and for two of them use meshgrid or ndgrid to create the appropriate matrices from the original vectors), and then solve for ‘s’ for each one (or each pair or combination) and then do a 2D or 3D plot of ‘s’ with respect to one or both of them.

In general, distance is a vector of a given magnitude and direction, so a negative distance would be a distance opposite the reference direction.

[Cesar Cardenas]

Thanks much for the explanation, I'm trying to get the result of s when z0 = 110, could show share it?

[Cesar Cardenas]

This equation:

can be typed in this way?

format long g
nz = @(z) n(z0)*(Tz0/Tz)^(1+g).*exp(-uc(z))    

[Cesar Cardenas]

just for clarification, when I run the code to solve for "s" I get:

Eqn =

s == -log((T_inf - Tz)/(T_inf - T_z0))/(z - z_0)

Eqn =

s == -(1.0*log((T_inf - 1.0*Tz)/(T_inf - 1.0*T_z0)))/(z - 1.0*z_0)

s =

-(1.0*log((T_inf - 1.0*Tz)/(T_inf - 1.0*T_z0)))/(z - 1.0*z_0)

but I do not get these s values

[Star Strider]

With respect to this Comment, that is a recursive call to ‘ni()’ so no. I have no idea how to code that. I am also not certain how to evaluate (that appears to be a constant), since it requires the recursive function as well.

With respect to this Comment, that is correct for that parameter set, and the same result I got for it.

[Cesar Cardenas]

Thank you, however when I run the code I can not get the numerial values, I just see the equations as shown in this Comment, so how could I obtain or see the numbers that are computed? Thanks

[Star Strider]

As always, my pleasure!

I am not certain what you are referring to. If you want to convert them to double values, use the double function.

[Cesar Cardenas]

now I understand better how the code works and I can finally see the value of "s" on my workspace. Thanks much.

[Cesar Cardenas]

Regarding this Comment , the derivate values is -236.40788......right?? which instruction calculate it?

D1Tz = subs(D1Tz) or D1Tz = vpa(D1Tz)?? or what is the difference between them??

Thank you

[Star Strider]

Correct.

The first assignment (with the subs call) creates it. In an earlier version of the code the first one originally produced a rational fraction, so the vpa call was necessasry to prodice a decimal fraction result. They are both correct, and the only difference may be the interim result, with the first returning a rational fraction and the second returning a decimal fraction. They are numerically the same.

[Cesar Cardenas]

Right thank you. If I wanted to find dT/dz using a value below z0=110 or above, the code would have to be adjusted or it can be done just by changing z0?? not sure if that makes sense. Thank you.

[Star Strider]

As always, my pleasure!

It would be necessary to do essentially everything in my earlier Comment with the appropriate parameters. If the parameters are different, just change them and run that code with them.

A more straightforward way would be to use ‘sfcn’ and ‘D1Tzfcn’ that I have introduced at the end of this Comment. They are anonymous functons that do not require the Symbolic Math Toolbox and can be run from base MATLAB. All the arguments of course must be present in the workspace when the functions are called.

syms s Tz T_inf T_z0 z z_0
Eqn = Tz == T_inf - (T_inf - T_z0) * exp(-s*(z-z_0))

Eqn = 

D1Tz = diff(rhs(Eqn),z)

D1Tz = 

Eqn = isolate(Eqn, s)

Eqn = 

% Tz = 200;
% z_0 = 120;
% T_inf = 651;
% T_z0 = 350;
% z = 120.7714;
% Eqn = subs(Eqn)
Eqn = vpa(Eqn)

Eqn = 

s = rhs(Eqn)

s = 

D1Tz = subs(D1Tz)

D1Tz = 

D1Tz = vpa(D1Tz)

D1Tz = 

sfcn = matlabFunction(rhs(Eqn))                         % Calculates 's'

sfcn = function_handle with value:
    @(T_inf,T_z0,Tz,z,z_0)(log((T_inf-Tz.*1.0)./(T_inf-T_z0.*1.0)).*-1.0)./(z-z_0.*1.0)

D1Tzfcn = matlabFunction(D1Tz)                          % Calculates 'dT(z)/dz'

D1Tzfcn = function_handle with value:
    @(T_inf,T_z0,Tz,z,z_0)(exp(log((T_inf-Tz.*1.0)./(T_inf-T_z0.*1.0)).*1.0).*log((T_inf-Tz.*1.0)./(T_inf-T_z0.*1.0)).*(T_inf-T_z0.*1.0).*-1.0)./(z-z_0.*1.0)

format longg
Tz = 200;
z_0 = 120;
T_inf = 651;
T_z0 = 350;
z = 120.7714;

s = sfcn(T_inf,T_z0,Tz,z,z_0)                           % Evaluate Function

s = 
        -0.524185992680584

dTdz = D1Tzfcn(T_inf,T_z0,Tz,z,z_0)                     % Evaluate Function

dTdz = 
         -236.407882698943

Using the anonymous functions to do the numeric calculations should make things easier. They are now independent of the Symbolic Math Toolbox. They should produce the same results with the same arguments, as demonostrated here, although with more limited precision (double instead of the symbolic extended precision, however that will likely not be a significant problem given their magnitudes).

.

[Cesar Cardenas]

Right, thanks so much for the clarification and explanation. Thank you.

[Star Strider]

As always, my pleasure!

[Cesar Cardenas]

Thanks much, just a quick question, based on the data in the image and this equation:

Tz = @(z) Tinf - (Tinf - Tz0).*exp(-s.*(z-z0));
Tz()

what would be a way to extend temperature upwards to find a temperature value at 200 km? Thanks.

[Star Strider]

I would not extrapolate those data. The atmosphere is not predictable enough to do that. Consider variations in the lapse rate, temperature inversions, and other potential nonlinearities. It is not possible to know, at least with any precision, what the atmosphere is doing outside the measured region.

I thought 112 km was in the ionosphere anyway. Is there even weather that high up?

[Cesar Cardenas]

I believe so, not sure though,, for this problem I would need to know the extend temperature values until 200 km for other calculations. Thanks.

[Star Strider]

You will have to measure them, not extrapolate them. I have no idea what ionospheric weather is.

The airplanes I flew never flew beyond the tropsphere, and the U.S. Federal Aviation Adminstrationi only requiered me to understand tropospheric weather.

[Cesar Cardenas]

Right thanks so much. Regarding this Comment, in that equation z is not included, just the change in altitude dz, I'm close to get the value I need but not still there, so, how could I adjust this code Comment?

Thanks for your help really appreciate it

[Cesar Cardenas]

just to clarify regarding my previous Comment, basically, I'm trying to find dT/dz like this:

would it be right?

Thanks much

T_1 = 184.829693;
T_2 = 160.199;
T_inf = 651;
T_z0 = 350;
z_1 = 109.94928741455;
z_2 = 106.9599609375;

s = (T_1 - T_2/z_1 - z_2)/(T_inf - T_z0)

[Star Strider]

That does not look correct to me.

[Walter Roberson]

looks like that would be the thermosphere https://en.m.wikipedia.org/wiki/Thermosphere

[Walter Roberson]

https://arxiv.org/pdf/1607.03370 seems relevant

[Star Strider]

I admit that I’ve never heard of the Thermosphere but then it doesn’t appear to affect tropospheric weather, so it was omitted from the aviation weather books I read. However, with a bit orf reading, some of this makes sense now. That context likely would have helped earlier.

Now I’m wondering if it affects long-distance shortwave (3 m to 30 m) radio communications (another interest of mine being amateur radio).

[Walter Roberson]

If I understood my earlier skimming properly, between 85 km and 115 km, the behaviour is dominated by mixing effects, so the temperature increase is largely unrelated to molecular weight, but above 115 km up to about 200 km, you get more and more stratification by molecular weight, and so you need a composite model that mixes between the 85/115 behaviour (that you have the data on) and other behaviour (which you do not have the data on.) The Bates-Walker model, if I understood correctly, is then considered valid between 85 and 115 but has to be ramped down above 115 km. And that the Bates-Walker model is old enough that those other effects were not known about, and the newer effects are still being studied so the newer models are still somewhat ad-hoc. If I understood the reading correctly, Bates-Walker is still used above 115 km by some people, but people are getting less and less satisfied with it.

All in all, it seems questionable to use the 88/113 data in the table to extrapolate up to 200 -- something that would have been acceptable in past times, and might be acceptable for a homework assignment these days, but not considered good enough for actual Science now ?

[Cesar Cardenas]

Yes this is basically a homework assignment in which the The Bates-Walker model is being applied. Thank you.

[Cesar Cardenas]

How can I code this equation with absolute value of dn/dz?

Hn = abs((dndz)/n)^-1?? not sure

Thanks much.

[Walter Roberson]

? The first one is a partial derivative, which could be positive or negative; the first one does not take magnitude in any way.

The second one might be a partial derivative or might be a simple derivative, and does appear to take magnitude.

If you are trying to implement the first form then you should not use abs()

[Cesar Cardenas]

Thanks. I think it is a simple derivate,,the first one maybe is not right in the notes I have.

[Walter Roberson]

Is it possible that n might be non-scalar? If so then the two forms could have very different meaning.

[Cesar Cardenas]

Thanks much. I'm trying to establish if it is a partial o simple derivative, it seems for sure it is a simple one.

Additionally, how could code and solve for qmax this equation:?? Thanks much.

[Walter Roberson]

If you know everything else then it looks like it would just be a division.

[Cesar Cardenas]

Yes, Thanks, so, it would be like this:?

q_max = q(h)/exp*(1 - (h - h_max)/H - sec(x)*exp(-(h-h_max)/H)

Thank you

[Walter Roberson]

Yes. Just remember to vectorize as appropriate.

[Cesar Cardenas]

Right thanks so much, as you can see in the previous Comment, q is q(h), (h is altitude) so how can I include it as a function of h in this equation, I have trying but I get errors. Thanks

q = 3.5;
h = 150;
h_max = 500;
H = 35.6757;
x = 15;

%format long g
q_max = q(h) ./exp((1 - (h - h_max)./H - sec(x)).*exp(-(h-h_max)./H))
q(150)
 

[Walter Roberson]

You don't.

I said earlier "If you know everything else", but you are saying now that you do not know everything else. You are trying to define in terms of and you are trying to define in terms of .

[Cesar Cardenas]

Right thank you. Maybe I was not so clear. For the problem I'm trying to solve, it is required to solve q(h) in terms of qmax and then the other way, based on the q(h) value obtained to solve for qmax, just not sure it is right?

[Walter Roberson]

Okay, the code you use in https://www.mathworks.com/matlabcentral/answers/1647255-how-can-i-write-and-solve-this-equation-as-a-function#comment_1994420 what errors do you encounter?

Perhaps you should be defining

q_max = @(h) q(h) ./exp((1 - (h - h_max)./H - sec(x)).*exp(-(h-h_max)./H))

[Cesar Cardenas]

Thanks, basically, for q(h) I have this:

I do not have errors and it runs well when the other parameters are defined.

q_max = 0.5;
%h = 139.6;
h_max = 198.369;
H = 35.6757;
x = 15;

format long g
qh = @(h) q_max.* exp((1 - (h - h_max)./H - sec(x)).*exp(-(h-h_max)./H));
qh (198)

I have this other equation for q(h) as in this case it depends on h and x, at first I get some erros but now it seems to run well, but sure it would be well defined?

q_max = 0.5;
%h = 198;
h_max = 198.369;
H = 35.6757;
%x = 15;

format long g
qh = @(h,x) q_max.*exp((1 - (h - h_max)./H - sec(x)).*exp(-(h-h_max)./H));
qh (198,15)

and for q_max as you mention in this comment I defined like this for instance for h = 100 and get errors

format long g
q_max = @(h) q(h) ./exp((1 - (h - h_max)./H - sec(x)).*exp(-(h-h_max)./H))
q_max(100)

Index exceeds the number of array elements. Index must not exceed 1.

Thanks much.

[Walter Roberson]

q is a variable that is an array; it is not a function at that point.

We do not see what q is in your code.

qh = @(h,x) q_max.*exp((1 - (h - h_max)./H - sec(x)).*exp(-(h-h_max)./H));

Okay, you made qh a function of h and x. That potentially makes sense

q_max = @(h) q(h) ./exp((1 - (h - h_max)./H - sec(x)).*exp(-(h-h_max)./H))

You made q_max a function of h but not of x . That feels a bit inconsistent after the definition for qh

Answer 2

[Cesar Cardenas]

Not sure what I'm not doing right as I'm not getting any change when the value is varied?? any help will be appreciated. Thanks.

net = 9.37e4;
T = 183;
mi= 5.31e-26;
k = 1.38e-23;
cs = 1.1304e-18;


format long g
Kz = @(z) cs*net*sqrt(8*k*T)/3.14*mi;
Kz(433)

[Walter Roberson]

Kz = @(z) cs*net*sqrt(8*k*T)/3.14*mi;

That function accepts 0 or 1 parameters, and ignores the parameter and computes a constant.

Note: we recommend you use pi instead of 3.14 unless you have particular reason to approximate the value.

Answer 3

[Cesar Cardenas]

Right thanks, how can I change it to accept any value? Thanks much

[Star Strider]

The function argument is ‘z’

Kz = @(z) cs*net*sqrt(8*k*T)/3.14*mi;

however ‘z’ nowhere appears in the expression to be evaluated.

What do you want it to be a function of?

[Cesar Cardenas]

Right, thanks, as I understand the problem to be solved, it would have to be a function of z, since the idea is to find the value of K at different z values, (altitudes). Thanks

[Walter Roberson]

Does cs vary with z? Does net vary with z? Does k vary with z? Does T vary with z? Does mi vary with z? Does the precision of pi that you use vary with z?

[Cesar Cardenas]

Thanks. Kz is proportional to this equation: (as seen in the image), so the idea is to find a z value where Kz maximizes. So, I would think (not sure) that v1,2 changes with altitude. Thanks.

[Star Strider]

What is ?

[Walter Roberson]

I do not not see any z in the ν definition. ν might change with altitude but which terms are constant with altitude and which vary with altitude?

[Cesar Cardenas]

As always thanks. Now, I think I have a better understanding of the problem. Actually Kz itself does not have to be calculated, just v1,2. What changes with altitude is n. The idea is to find the max. value of v1,2 as n varies within the 105-109 km range. I already coded v1,2, but not sure how I can complete/arrange the code to find a maximum value with 105-109. As always any help will be appreciated. Thanks.

cs = 1.1304e-18;
n = 9.37e4;
T = 183;
k = 1.38e-23;
mi= 5.31e-26;

v12 = cs*n*sqrt(8*k*T)/pi*mi

[Star Strider]

‘What changes with altitude is n. The idea is to find the max. value of v1,2 as n varies within the 105-109 km range.’

The ‘v12’ variable is apparently a linear function of ‘n’ so the maximum value of ‘v12’ will be at the maximum value of ‘n’, unless other variables in the equation also change with altitude.

[Cesar Cardenas]

Thanks much. I already solved this problem. What would be a nice approach for this one: (not sure about it).

As always any help will be greatly appreciated. Thanks.

[Walter Roberson]

That appears to be an unrelated question.

[Cesar Cardenas]

Thanks much. I'll check it more carefully. Aside this, wondering if this equation would be well coded? As always, thanks much.

G = 74.7500e9;
b = 0.445e-9;
v = 0.35;
d = 10e10;

Ur = G*b^2./(10) + G*b^2./(4*pi*(1-v))*log(d^(-1./2)./(5*b))

[Star Strider]

I don’t believe it is.

There’s a term in the denominator of the second term, that is then mulitplied by so they should cancel out. Either that, or there is a typographical error in the second term.

[Walter Roberson]

I was curious about which power operation would be fastest. The three possibilities are ^(-1./2) (matrix power), .^(-1./2) (power), and 1/sqrt() (potentially using hardware square root).

Conclusions are a bit tricky.

If you look at the below code, notice I discard the first few samples. When you do not do that, then right at the beginning you get some timings that are much higher.

If you look at the graph below, notice there is an isolated large peak for ^ and another for .^ . In my tests up to N = 1000, it was not unusual to see those at isolated places, especially for ^ -- for some reason the ^ operator seems more "bursty" than the others, which I have no explanation for at the moment (because it should devolve to the same code as .^ for the case of scalars). If you use movmedian() with a mean of 3 then those peaks disappear.

You can see from the bottom half of the plot that the timings are pretty variable for the three different operations, with none of them being a clear winner. So we have to look at something like the mean() to talk about the statistical behaviour.

In my tests, most of the time ^ was slightly slower on average. That is not unexpected, as it has to do an extra step to determine whether to use scalar behaviour. Occasionally .^ shows up slower on average, which I suspect is due to those occasional large bursts.

In my tests, 1/sqrt() was almost always faster on average; with N = 1000, this becomes more obvious.

%initialize
N = 300;
t_mpower = zeros(N,1); t_power = zeros(N,1); t_sqrt = zeros(N,1);

%time
for K = 1 : N; d = rand(); start = tic; d^(-1./2); t = toc(start); t_mpower(K) = t; end
for K = 1 : N; d = rand(); start = tic; d.^(-1./2); t = toc(start); t_power(K) = t; end
for K = 1 : N; d = rand(); start = tic; 1./sqrt(d); t = toc(start); t_sqrt(K) = t; end

%discard "warm-up"
t_mpower(1:9) = []; t_power(1:9) = []; t_sqrt(1:9) = [];

%display
format long g
[mean(t_mpower), mean(t_power), mean(t_sqrt)]

ans = 1×3

1.0e+00 * 6.04810996563573e-07 5.60137457044673e-07 4.7766323024055e-07

plot(([t_mpower, t_power, t_sqrt]))
legend({'\^', '.\^', '1/sqrt'})

[Walter Roberson]

log(d^(-1./2)./(5*b))

Note that that should be the same as

-log(5 .* b .* sqrt(d))

which could potentially be further decomposed (might be able to pre-calculate parts of that, depending what is varying.)

[Cesar Cardenas]

Thanks much, I missed that. would it be right now? Thanks much.

G = 74.7500e9;
b = 0.445e-9;
v = 0.35;
d = 10e10;

Ur = G*b^2./(10) + G*b^2./(4*pi)*(log(d^(-1./2)./(5*b)))

[Cesar Cardenas]

Thanks much interesting approach.

[Cesar Cardenas]

Thank you, would there be any way to relate the below code for a possible solution to this previous Comment?? As always any help will be greatly appreciated. Thanks

close all
clear all
clc
x=linspace(-5*10^9,5*10^9,101);
y=linspace(-5*10^9,-1*10^9,101);
G=27.7*10^9;
b=0.286*10^-9;
v=0.334;
[x,y]=meshgrid(x,y);
sxx=(-G*b)/(2*pi*(1-v))*y.*(3*x.^2+y.^2)./(x.^2+y.^2).^2;
syy=(G*b)/(2*pi*(1-v))*y.*(x.^2-y.^2)./(x.^2+y.^2).^2;
sxy=(G*b)/(2*pi*(1-v))*x.*(x.^2-y.^2)./(x.^2+y.^2).^2;
surf(x,y,sxx);
title('Plot of sxx');
figure;
surf(x,y,syy);
title('Plot of syy');
figure;
surf(x,y,sxy);
title('Plot of sxy');

[Star Strider]

It appears to be coded correctly (although not the way I would code it). That aside, the plot appear to be appropriate, although stress fields not an area of my expertise.

x=linspace(-5*10^9,5*10^9,101);
y=linspace(-5*10^9,-1*10^9,101);
G=27.7*10^9;
b=0.286*10^-9;
v=0.334;
[x,y]=meshgrid(x,y);
sxx=(-G*b)/(2*pi*(1-v))*y.*(3*x.^2+y.^2)./(x.^2+y.^2).^2;
syy=(G*b)/(2*pi*(1-v))*y.*(x.^2-y.^2)./(x.^2+y.^2).^2;
sxy=(G*b)/(2*pi*(1-v))*x.*(x.^2-y.^2)./(x.^2+y.^2).^2;

figure
surfc(x,y,sxx);
title('Plot of sxx');

figure;
surfc(x,y,syy);
title('Plot of syy');

figure;
surfc(x,y,sxy);
title('Plot of sxy');

.

[Cesar Cardenas]

@Star Strider thanks much, how it can be plotted by using the filled contour plot options? Thanks much.

[Star Strider]

As always, my pleasure!

Replace the surf (or surfc) calls with contourf calls. The arguments remain the same.

x=linspace(-5*10^9,5*10^9,101);
y=linspace(-5*10^9,-1*10^9,101);
G=27.7*10^9;
b=0.286*10^-9;
v=0.334;
[x,y]=meshgrid(x,y);
sxx=(-G*b)/(2*pi*(1-v))*y.*(3*x.^2+y.^2)./(x.^2+y.^2).^2;
syy=(G*b)/(2*pi*(1-v))*y.*(x.^2-y.^2)./(x.^2+y.^2).^2;
sxy=(G*b)/(2*pi*(1-v))*x.*(x.^2-y.^2)./(x.^2+y.^2).^2;

figure
contourf(x,y,sxx);
title('Plot of sxx');

figure;
contourf(x,y,syy);
title('Plot of syy');

figure;
contourf(x,y,sxy);
title('Plot of sxy');

.

[Cesar Cardenas]

@Star Strider Thanks so much, how the coordinate range [-50,50] and contour level 10 can be set up here?? Thanks much.

[Star Strider]

I don’t understand.

My best guess —

x=linspace(-5*10^9,5*10^9,101);
y=linspace(-5*10^9,-1*10^9,101);
G=27.7*10^9;
b=0.286*10^-9;
v=0.334;
[x,y]=meshgrid(x,y);
sxx=(-G*b)/(2*pi*(1-v))*y.*(3*x.^2+y.^2)./(x.^2+y.^2).^2;
syy=(G*b)/(2*pi*(1-v))*y.*(x.^2-y.^2)./(x.^2+y.^2).^2;
sxy=(G*b)/(2*pi*(1-v))*x.*(x.^2-y.^2)./(x.^2+y.^2).^2;

figure
contourf(x,y,sxx, 10);
xlim([-50  50])
colormap(turbo(10))
title('Plot of sxx');

figure;
contourf(x,y,syy, 10);
xlim([-50  50])
colormap(turbo(10))
title('Plot of syy');

figure;
contourf(x,y,sxy, 10);
xlim([-50  50])
colormap(turbo(10))
title('Plot of sxy');

.

[Cesar Cardenas]

Thanks much. Maybe I was not so clear, the idea is to use the contour plot in the coordinate range of [-50, 50] and I think the solution would be similar to the one in the previous comment. Not sure how to do it.

Thanks so much

[Star Strider]

I do not understand ‘coordinate range of [-50, 50] ’

What coordinates does this refer to?

Answer 4

[Cesar Cardenas]

Thanks much I think your approach here comment. is pretty similar to the below image. I think the coordinate range may refer to something like this graph...not sure though. Thanks much.

[Star Strider]

I’m still lost.

This is not an area of my expertise.

Perhaps re-defining ‘x’ and ‘y’ in the earlier code will do what you want.

Answer 5

[Cesar Cardenas]

How can I convert this value 0.6667 to a fraction? Thanks

[Star Strider]

rats!

d = 0.6667;
strlen = 3;
f = rats(d,strlen)

f = ' 2/3'

.

[Cesar Cardenas]

Thanks much, and for instance if I want to convert a fraction to a radical? , eg, 1/2 to radical? Than you.

[Star Strider]

Does ‘radical’ mean ‘root’ or decimal fraction?

[Cesar Cardenas]

root, thanks much

[Walter Roberson]

sym() applied to a floating point number attempts to find fractions and square roots and pi

[Star Strider]

@Walter Roberson — Thank you! (I was off editing another answer.)

@Cesar Cardenas —

Still not certain.

The root of is

format long
sqrt(1/2)

ans = 
   0.707106781186548

The square of is of course .

.

[Cesar Cardenas]

Thanks much. That wil help. Additionally, I'm trying to code this:

This is what I'm getting. not sure how I can get a single value for rl?? not this vector..thanks much. As always any help will be greatly appreciated.

m = 9.11e-31;
q = -1.6e-19;
v = [400 200 300]; %velocity
B = [2 30 0];%uniform magnetic field

FB = q*cross(v,B*10^-3)

rl = m.*v./abs(q).*FB1

rl =

   1.0e-26 *

    0.3280   -0.0109   -0.3170

[Star Strider]

I have no idea.

This is far from the original topic, and not an area of my expertise.

[Walter Roberson]

@Cesar Cardenas

v is a vector, so m.*v is a vector and v./abs(q) is a vector. FB1 is a vector and vector .* vector is a vector.

not sure how I can get a single value for rl??

Look at your equation. The F is underlined, indicating that the result is expected to be a vector, not a scalar.

[Cesar Cardenas]

Right thanks much, I'll check it carefully again as I think it should be a scalar. Additionally, I'm trying to code this equation:

This is my initial approach but since it is a function of z, not sure how to arrange it?? as always any help will be appreciated. Thank you.

lsr = 3e10; 
ne = 1.3e10;
ni = 2e16:
z = 160;

Bz = lsr./ni.*ne 

[Cesar Cardenas]

How can I calculate the magnitude of this vector:

v = [2 30 0]*10^-3,, it is to 10^-3 not sure how to add it to the code

v = [2: 30: 0];
sv = v.* v;       %the vector with elements 
                  % as square of v's elements
dp = sum(sv);     % sum of squares -- the dot product
mag = sqrt(dp);   % magnitude
disp('Magnitude:'); 
disp(mag);

[Cesar Cardenas]

@Walter Roberson thank you regarding this comment. Basically here what I need to is find V perpendicular as shown here;

So first I did this:

%Magnitude of B
B = [2 30 0];%uniform magnetic field
norm(B*10^-3)

Then, this; to find V parallel,

dot(v,B)

So, would it be a nice approach? or how can I work this out? Thanks much