Solving/plotting a nonlinear equation for multiple values

Hello all,

I am a mechanical engineering student working on some Applied Thermodynamics homework and am having trouble getting the values I need to plot a nonlinear equation. I have (on paper) determined that my nonlinear equation is as follow:

2*sqrt(P)*((x)^1.5) + .10905*x- 0.0363 = 0

where P is a pressure and x happens to be a molecular concentration in this case. So, I can solve the above equation for set values of P using the following function file (when P = 1).

function [x,P] = myfun(x)
F = 2*sqrt(1)*((x)^1.5) + .10905*x- 0.03635;
end

and the commands:

x0 = 0.5;
[x,fval] = fsolve(@myfun,x0)

But what I want is to define P as an array like [0.1:0.1:100] and be able to solve all of those values at once instead of having to get each one individually. I was wondering if this was possible using the above approach or if I should be looking at another method. Thank you in advance for any guidance you can give me. Still kind of a Matlab novice if you couldn't already tell.

Answer 1

[Walter Roberson]

There are three solutions, two of which are only real-valued between P = 0 and P = 0.000132 (approximately). The third is real-valued for larger P as well.

The order of floating point calculations can make a big difference in functions such as these, sometimes changing the calculations completely. I have therefore converted the floating point values you gave into rational values; the below is in rational

((1/40000)*((581600000000*P^2 - 384240583*P + 29080000*(400000000*P - 528529)^(1/2) * P^(3/2))/P^(5/2))^(1/3) + (528529/40000) / (P*((581600000000*P^2 - 384240583*P + 29080000*(400000000*P - 528529)^(1/2)*P^(3/2)) / P^(5/2))^(1/3)) - (727/40000) / P^(1/2))^2

 (1/4652800000000)*727^(1/3)*((-1+I*3^(1/2))*P^(3/2)*727^(2/3) * ((800000000*P^2 - 528529*P + 40000*(400000000*P - 528529)^(1/2) * P^(3/2)) / P^(5/2))^(2/3) - 1454*P*727^(1/3)*((800000000*P^2 - 528529*P + 40000*(400000000*P - 528529)^(1/2)*P^(3/2)) / P^(5/2))^(1/3) - 528529*(1+I*3^(1/2))*P^(1/2))^2 / (((800000000*P^2 - 528529*P + 40000*(400000000*P - 528529)^(1/2)*P^(3/2))/P^(5/2))^(2/3)*P^3)

 (1/4652800000000)*((1+I*3^(1/2))*P^(3/2)*727^(2/3) * ((800000000*P^2 - 528529*P + 40000*(400000000*P - 528529)^(1/2)*P^(3/2))/P^(5/2))^(2/3) + 1454*P*727^(1/3)*((800000000*P^2 - 528529*P + 40000*(400000000*P - 528529)^(1/2)*P^(3/2)) / P^(5/2))^(1/3) - 528529*P^(1/2)*(-1+I*3^(1/2)))^2*727^(1/3) / (((800000000*P^2 - 528529*P + 40000*(400000000*P - 528529)^(1/2)*P^(3/2)) / P^(5/2))^(2/3)*P^3)

[Weston]

So, the third function you have listed there is equal to the x that was in the original equation? How did you get that solution and what is the I term? Sorry if these are rudimentary questions, I am just slightly confused. Thank you very much for the help.

[Walter Roberson]

Yes, all three of those are solutions to x from your original question. The second and third are valid up to about 0.000132; the first is valid for positive P.

I got the equations by converting the floating point values to rational, and then using a symbolic toolkit to solve() the equation.

[Weston]

Ah, I am not yet versed in the symbolic toolkit though it sounds like a lot more user friendly way to do this type of problem. Though I am still not sure what I is in the above solution.

[Walter Roberson]

Sorry, "I" is sqrt(-1)

You can essentially rewrite the equation into the form

a = .10905 / (2*sqrt(P))
b = -0.0363 / (2*sqrt(P))
x^(3/2)+a*x+b

and then solve() for x. The solution becomes fairly similar to the solution for cubics, but I do not myself understand how the solution is arrived at.

With a symbolic toolbox, you would go directly to solving x^(3/2)+a*x+b for x, and then subs() in the values for "a" and "b", then simplify(). Then you could matlabFunction() the result to get a vectorized function that takes in a vector of P and returns the corresponding x. You might need to filter out the complex values for x afterwards, depending how you want to handle that small blip of three solutions at low pressures.

[Weston]

Perfect. That makes a lot of sense. I have a much better idea of what's going on now. Thank you so much for your help.

Answer 2

[Alfonso Nieto-Castanon]

Perhaps something like:

 f = @(x,P)2*sqrt(P)*x^1.5 + .10905*x- 0.0363;
 P = .1:.1:100;
 x0 = 0.1;
 [x,fval] = arrayfun(@(P)fzero(@(x)f(x,P),x0), P);

 plot(P,x);