how to code equation

Hello, I'm trying to code this equation but I think it is not right. Not sure how to adjust the LHS to avoid errors. Any help will be greatly appreciated. Thanks.

psi - log(psi) = -3 * -4 * log(l/2) + 4 * log(E) + 2 * l/E 

[Walter Roberson]

See by the way https://www.mathworks.com/matlabcentral/answers/225896-how-can-i-calculate-ln-x-in-matlab-code#answer_1195970 in which I talk about how different computer languages represent natural logarithm in practice.

Answer 1

[Star Strider]

Using the Symbolic MAth Toolbox —

syms xi lambda psi  
Eqn = psi - log(psi) == -3 * -4 * log(lambda/2) + 4 * log(xi) + 2 * lambda/xi 

Eqn = 

You can then manipulate it symbolically and calculate with it symbolically.

.

[Cesar Cardenas]

Thank you for your reply. if psi = 1 and xi = lambda/2, how do I have to define them here? Thank you once again.

[Torsten]

syms xi lambda psi  
Eqn = psi - log(psi) == -3 * -4 * log(lambda/2) + 4 * log(xi) + 2 * lambda/xi ;
Eqn_subs = subs(Eqn,[psi,xi],[sym('1'),lambda/sym('2')])

Eqn_subs = 

lambda_num = solve(Eqn_subs,lambda)

ans = 

[Star Strider]

@Torsten — Thank you. (Away for a few hours here.)

[Cesar Cardenas]

Thank you so much for you help. For the values provided for psi = 1 and xi = lambda/2, just wondering how I can plot and get a curve like the ones shown here? I'm not a matlab expert and this is the first time I use the symbolic tool. Again, any help will be appreciated. Thank you.

[Star Strider]

The value of ξ seems to go from about 1 to about 90. Doing those substitutions as stated and solving for λ produces this —

syms xi lambda psi  
Eqn = psi - log(psi) == -3 * -4 * log(lambda/2) + 4 * log(xi) + 2 * lambda/xi;
Eqn = subs(Eqn, {xi,psi}, {lambda/2,1})

Eqn = 

Eqn = solve(Eqn, lambda)

Eqn = 

I am not following what you want to do, since this will not produce one curve, much less a family of curves, since it is not a function of any variable.

Please clarify.

.

[Walter Roberson]

The plot shows v(r) as a function of ξbut I do not see any v(r) in the formula ?

I am guessing that the 0.5*10^6 and so on are different values of ψ but you have been talking about ψ as 1 rather than as 1*10^6 ?

[Cesar Cardenas]

@Walter Roberson Yes, I think you're right based on what I've read. That's what I would need to plot. Any clue on how I could get it? Thanks much

[Star Strider]

What variable needs to vary and what is the range?

If more than one needs to vary, what are they and what are the ranges of each?

[Cesar Cardenas]

@Star Strider The idea would be to reproduce the plot shown in the previous image. So, the same ranges could be used. Not surre how to do it with the symbolic tool. As always, any help will be appreciated. Thanks

[Walter Roberson]

@Cesar Cardenas

None of the variables in psi - log(psi) = -3 * -4 * log(l/2) + 4 * log(E) + 2 * l/E are v(r) so we do not know what it is you want to plot.

We also need to know whether the values for ψ should be near 1 (as in your original question) or if they should be near 1 * 10^6 (as in your plot)

[Cesar Cardenas]

@Walter Roberson @Star Strider I have set up this alternative equation that would lead to a similar solution. The idea is to get a curve or possible solution as shown below. (not all of them). As stated before, with the symbolic option I'm not sure how I could get a plot like this? Any help will be appreciated. Thanks

R = 6.96e5;
r = 1e11;
u = 1.7e-2;
c = 40;
g % gravity

syms c u R r
Eqn = 1/2 * u^2 - c^2 * log(u) == 2*c^2 *log(r+g)* R^2/r

[Walter Roberson]

That equation is not even close for the specific input values you list.

In the symbolic form, I cannot tell what you intend to solve for.

Your equations do not involve V or or so we cannot tell what you are trying to do.

Q = @(V) sym(V);

R = Q(696)*Q(10)^3;
r = Q(1)*Q(10)^11;
u = Q(17)*Q(10)^(-3);
c = Q(40);
g = Q(981)*Q(10)^(-2); % gravity

%syms c u R r
Eqn = 1/2 * u^2 - c^2 * log(u) == 2*c^2 *log(r+g)* R^2/r

Eqn = 

lhs(Eqn) - rhs(Eqn)

ans = 

vpa(ans)

ans = 

[Cesar Cardenas]

Thank you for your reply and help. I looked into this a bit more. I set up another similar equation. What I need is to plot just or at least one curve as the shown below. Not sure how to plug the uc and rc values in Eqn and plot it. As always any help will be appreciated. Thanks much.

k = physconst('Boltzmann');
T = 1.2e6;
mH = 1.66e-27;
G = 6.67e-11;
Ms = 1.99e30;

uc = sqrt(2 * k * T/mH);
rc = G * Ms * mH/4 * k * T;

syms uc u r rc
Eqn = (u/uc)^2 - 2 * log(u/uc) == 4 * log(r/rc) + 4 * log(rc/r) - 3 

[Walter Roberson]

uc = sqrt(2 * k * T/mH);
rc = G * Ms * mH/4 * k * T;

Those are constants.

syms uc u r rc

That has the effect of

uc = sym('uc');
rc = sym('rc');

which would overwrite the numeric scalars with symbolic variable names.

You have a single equation in 3 variables. You can solve for any one of them, but you would not be able to get a 2D plot out of it.

[Cesar Cardenas]

Right, thank you, is there any other option to do this? Thanks

[Star Strider]

The arguments you want to plot agiainst must be vectors, since at least one variable must be a vector in order to plot the evaluation of the function against it.

Probably the easiest way to create vectors that span two limits (a minimum value to a maximum value) is to use the linspace function. It is then also straightforward to be certain all of them have the same lengths, regardless of the range of values that they span. The functional plotting functions such as fplot, fsurf, and the others can provide one or more of those vectors themselves to produce a plot.

[Cesar Cardenas]

Right thank you. Could you give at least any clue on how to set up the code for this? I'm not a matlab expert. Thank you.

[Walter Roberson]

Q = @(v) sym(v);

k = Q(physconst('Boltzmann'));
T = Q(1.2)*(1e6);
mH = Q(1.66)*(1e-27);
G = Q(6.67)*(1e-11);
Ms = Q(1.99)*(1e30);

uc = sqrt(2 * k * T/mH);
rc = G * Ms * mH/4 * k * T;

syms u r real
Eqn = (u/uc)^2 - 2 * log(u/uc) == 4 * log(r/rc) + 4 * log(rc/r) - 3 

Eqn = 

simplify(Eqn)

ans = 

sol = simplify(solve(Eqn, u))

Warning: Possibly spurious solutions.

Warning: Solutions are only valid under certain conditions. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

sol = 

sol1 = sol(1);
sol1r = real(sol1);
sol1i = imag(sol1);
tiledlayout('flow');
nexttile(); fplot(sol1r, [0.01 1]); title('real component')
nexttile(); fplot(sol1i, [0.01 1]); title('imaginary component')

Look at the right hand side of your Eqn. You have

4 * log(r/rc) + 4 * log(rc/r) - 3

When r and rc are positive and non-zero then log(r/rc) = -log(rc/r) and the log terms on the right hand side cancel.

Eqn1 = (u/uc)^2 - 2 * log(u/uc) == - 3 

Eqn1 = 

solve(Eqn1, u)

 
ans =
 
Empty sym: 0-by-1
 

syms x
solve( x^2 - 2*log(x) == -3 )

ans = 

vpa(ans)

ans = 

so u/uc would have to be complex-valued for there to be a solution.