**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

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

20 views (last 30 days)

Show older comments

##### 1 Comment

Benjamin Thompson
on 10 Feb 2022

### Accepted Answer

Star Strider
on 10 Feb 2022

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)

.

##### 53 Comments

Cesar Cardenas
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

dTdz = @(z) appropriate calculation of derivative

s = dTdz(z0) / (Tinf - Tz0)

Star Strider
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

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
on 11 Feb 2022

Star Strider
on 11 Feb 2022

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
on 11 Feb 2022

Cesar Cardenas
on 11 Feb 2022

This equation:

can be typed in this way?

format long g

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

Cesar Cardenas
on 12 Feb 2022

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
on 12 Feb 2022

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
on 12 Feb 2022

Edited: Cesar Cardenas
on 12 Feb 2022

Star Strider
on 12 Feb 2022

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
on 12 Feb 2022

Cesar Cardenas
on 12 Feb 2022

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
on 12 Feb 2022

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
on 12 Feb 2022

Star Strider
on 12 Feb 2022

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
on 14 Feb 2022

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

Cesar Cardenas
on 15 Feb 2022

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
on 15 Feb 2022

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
on 15 Feb 2022

Star Strider
on 15 Feb 2022

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
on 15 Feb 2022

Cesar Cardenas
on 15 Feb 2022

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)

Walter Roberson
on 16 Feb 2022

Edited: Walter Roberson
on 16 Feb 2022

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

Star Strider
on 16 Feb 2022

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
on 16 Feb 2022

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
on 16 Feb 2022

Cesar Cardenas
on 17 Feb 2022

Edited: Cesar Cardenas
on 17 Feb 2022

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

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

Thanks much.

Walter Roberson
on 17 Feb 2022

? 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
on 18 Feb 2022

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

Walter Roberson
on 18 Feb 2022

Cesar Cardenas
on 18 Feb 2022

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
on 18 Feb 2022

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

Cesar Cardenas
on 18 Feb 2022

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

Cesar Cardenas
on 19 Feb 2022

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
on 19 Feb 2022

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
on 19 Feb 2022

Edited: Cesar Cardenas
on 19 Feb 2022

Walter Roberson
on 20 Feb 2022

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
on 20 Feb 2022

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
on 20 Feb 2022

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

### More Answers (4)

Cesar Cardenas
on 16 Mar 2022

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)

##### 1 Comment

Walter Roberson
on 16 Mar 2022

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.

Cesar Cardenas
on 16 Mar 2022

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

##### 24 Comments

Star Strider
on 16 Mar 2022

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
on 16 Mar 2022

Walter Roberson
on 16 Mar 2022

Cesar Cardenas
on 16 Mar 2022

Walter Roberson
on 17 Mar 2022

Cesar Cardenas
on 17 Mar 2022

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
on 17 Mar 2022

‘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
on 18 Mar 2022

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.

Cesar Cardenas
on 20 Mar 2022

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
on 20 Mar 2022

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
on 20 Mar 2022

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
on 20 Mar 2022

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
on 20 Mar 2022

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
on 20 Mar 2022

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
on 20 Mar 2022

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
on 21 Mar 2022

Star Strider
on 21 Mar 2022

As always, my pleasure!

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
on 21 Mar 2022

Star Strider
on 21 Mar 2022

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
on 21 Mar 2022

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
on 21 Mar 2022

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

What coordinates does this refer to?

Cesar Cardenas
on 22 Mar 2022

##### 1 Comment

Star Strider
on 22 Mar 2022

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.

Cesar Cardenas
on 27 Mar 2022

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

##### 12 Comments

Star Strider
on 27 Mar 2022

Cesar Cardenas
on 27 Mar 2022

Walter Roberson
on 27 Mar 2022

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

Star Strider
on 27 Mar 2022

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

Still not certain.

The root of is

format long

sqrt(1/2)

ans =

0.707106781186548

The square of is of course .

.

Cesar Cardenas
on 29 Mar 2022

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
on 29 Mar 2022

I have no idea.

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

Walter Roberson
on 29 Mar 2022

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
on 31 Mar 2022

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
on 31 Mar 2022

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
on 31 Mar 2022

@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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)