I have a nonlinear equation with a symbolic variable and cant solve it.

Hi, I have to plot a phase diagram of an implicit equation, but in there lies a function p that depends nonlinearly on one of the variables (alpha), so I am trying to solve the nonlinear equation in terms of my variable alpha. However, I only get a solution in terms of a z2 variable. Can anybody help? Here is the code

n = 2.1;
alpha0 = 0;
syms p alpha
eqn = p^(n+1) - alpha0*p^n + p -(alpha0 + alpha) == 0;
sol = solve(eqn,p);

Warning: Solutions are parameterized by the symbols: z2. To include parameters and conditions in the solution, specify the 'ReturnConditions' value as 'true'.

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

For n = 2 works just fine, but not for n = 2.1, which prompts my problem.

Answer 1

[John D'Errico]

Look carefully at your problem.

n = 2.1;
alpha0 = 0;
syms p alpha
eqn = p^(n+1) - alpha0*p^n + p -(alpha0 + alpha) == 0

eqn = 

So you have sort of, almost kind of, a cubic polynomial, but not really. Raising numbers to non-integrer powers causes problems. There is NO analytical solution. Effectively, we can think of this as a very high order polynomial problem, of degree 31 in fact.) And polynomials higher than order 4 have no solution in algebraic form. This has been known for well over 100 years.

Is there something you could do? If you have the value of alpha, then you could simple use fzero. For example, you might do this:

Psol = @(alpha) fzero(@(p) p^(n+1) - alpha0*p^n + p -(alpha0 + alpha),1)

Psol = function_handle with value:
    @(alpha)fzero(@(p)p^(n+1)-alpha0*p^n+p-(alpha0+alpha),1)

Psol(3)

ans = 1.2072

So plug in any value of alpha, and fzero will (hopefully) return a solution. You could also use vpasolve, which will survive as long as alpha is given in numerical form. Or fsolve, etc. Another way of looking at this is in the form of a plot. (If we recognize that whenever p is less than zero, all sorts of hell will break lose in terms of complex variables, then we can see that p should be positive. So write the problem by isolating alpha from the rest.

alpha = p + p^3.1

Now we can plot alpha, as a function of p.

fplot(@(p) p + p.^(3.1),[0,3])
xlabel p
ylabel alpha
grid on

Given any value of alpha on the vertical axis, you wish to solve for p. Again, there is NO algebraic solution to this. Yes, if you wanted, you could approximate the inverse relationship. What you want to see is that same relation, but with the axes flipped. And you could easily enough do that. For example...

pint = linspace(0,10);
alphaint = pint + pint.^3.1;
plot(alphaint,pint,'-')
xlabel alpha
ylabel p
grid on

You could now easily enough fit those points with a spline. For example...

pspl = spline(alphaint,pint);
fnval(pspl,5)

ans = 1.4982

So when alpha is 5, p would now directly be predicted as 1.4982. Of course, this is just an interpolation. And again, there is no algebraic function you can write down. Sorry.

[Álvaro]

The span of alpha for the figure I want to plot is [1 10^4], so there exists the posibility of running the equation for all possible values of alpha, let's say for x = linspace(1,10^4,10^5). However, when running this in a "for" loop, it almost instantaneously goes crazy as the initial condition is probably off (I think, though it says it is due to encountering complex solutions). The code is like this

n = 2.1;
alpha0 = 0;
x = linspace(1,10^4,10^5);
A = zeros(10^5,1);
Psol = @(alpha) fzero(@(p) p^(n+1) - alpha0*p^n + p -(alpha0 + alpha),1);


for i = 1:length(x)
    A(i,1) = Psol(i);
end

Incorrect number or types of inputs or outputs for function solve.

Error in solution>@(alpha)solve(@(p)p^(n+1)-alpha0*p^n+p-(alpha0+alpha),1) (line 5)
Psol = @(alpha) solve(@(p) p^(n+1) - alpha0*p^n + p -(alpha0 + alpha),1);

[Álvaro]

I have also tried with the vpasolve, which has finally worked with the following code

n = 2.1;
alpha0 = 0;
x = linspace(1,10^4,10^5);
A = zeros(10^5,1);
syms p
for i = 1:length(x)
    eq = p^(n+1) - alpha0*p^n + p -(alpha0 + i) == 0;
    A(i,1) = vpasolve(eq,p);
end

My computer has been working for 30 minutes but I finally have an array for all the values of p from alpha = 1 to alpha = 10^4. I do not know whether the values are correct or how to implement this result onto a fimplicit equation, but I will manage. Thanks a lot!

[John D'Errico]

Honestly, you are going the wrong way. First, how far out does p go, if alpha is 1e4?

n = 2.1;
alpha0 = 0;
x = linspace(1,10^4,10^5);
A = zeros(10^5,1);

Note that I used fzero slightly different here, giving it a bracket. As long as the upper end is high enough, it will have no problem.

Psol = @(alpha) fzero(@(p) p^(n+1) - alpha0*p^n + p -(alpha0 + alpha),[0,100]);
Psol(1e4)

ans = 19.5007

Ok, so we want to go as high as 19.5007, or so.

pint = linspace(0,20,10000); % a nice fine spacing
alphaint = pint.^(n+1) - alpha0*pint.^n + pint -alpha0;

Now we have a set of values that relate alpha and p.

plot(alphaint,pint,'.')
xlabel alpha
ylabel p
grid on

And the nice thing is, this was virtually instantaneous. Now you can use any interpolant to recover p, as a function of alpha.

Perhaps even better yet, see this relation is fairly well-behaved if we plot it in loglog form

loglog(alphaint,pint,'-')
xlabel alpha
ylabel p
grid on

And that suggests you have a simple scheme that will work very nicely. Interpolate log(p), as a function of log(alpha). Even a linear interpolant will sufficent there. Then exponentiate the result. something like this:

ppred = @(alpha) exp(interp1(log(alphaint(2:end)),log(pint(2:end)),log(alpha)));

ppred(3000)

ans = 13.2140

Note that I dropped the first point in those vectors, since it will be zero, and log(0) is a problem.