How to find the perturbation of $x^2 − 1 = \epsilon x$?

Question

Is there a function in Mathematica that can be used to find the perturbation solution of an equation like $x^2 − 1 = \epsilon \,x$, $x − 2 = \epsilon \cosh(x)$ or $x^2 − 1 = \epsilon\, e^x$?

Answer 1

Decide up to which power you would like to expand:

pow = 4;

Let's do one of the equations you mentioned as an example (bring all terms to one side and save as a single expression, eq in this case):

eq = x - 2 - e Cosh[x];

Write an ansatz for the x solution with unknown coefficients:

ansatz = Sum[a[i] e^i, {i, 0, pow}];

Substitute the ansatz into your expression and do a series expansion:

expand = Series[eq /. x -> ansatz, {e, 0, pow}];

Solve the constraints of overall factors vanishing:

vars = Table[a[i], {i, 0, pow}];
sols = Solve[expand == 0, vars][[1]] // Simplify // Quiet;

Finally, insert the solution into your ansatz to obtain the result:

res = ansatz /. sols

2 + e Cosh[2] + e^2 Cosh[2] Sqrt[-1 + Cosh[2]^2] + 1/3 e^4 Cosh[2] Sqrt[-1 + Cosh[2]^2] (-3 + 8 Cosh[2]^2) + e^3 (-Cosh[2] + (3 Cosh[2]^3)/2)

Don't forget to test numerically, whether your result is actually correct at the end:

enum = 10^-10;
xnum = x /.FindRoot[(eq /. e -> enum) == 0, {x, 2}, WorkingPrecision -> 60];
(res /. e -> enum) - xnum

-3.627605747396*10^-47

This shows that expanding to fourth order with an e= 10^-10 indeed consistently matches the result up to about 10^-50 accuracy, so the expansion was correct. Rinse and repeat for the other examples.

PS: In cases where your equation admits several solutions you might have to be a little bit more careful, but the principle still stays the same.

Answer 2

For these simple examples you could use InverseSeries. For example:

InverseSeries[Series[(x^2 - 1)/x, {x, 1, 10}], e]
InverseSeries[Series[(x - 2)/Cosh[x], {x, 2, 10}], e]
InverseSeries[Series[(x^2 - 1)/E^x, {x, 1, 10}], e]

You need to solve for e, and then do a series around the value of x when e is zero.

Answer 3

The new in M12 function AsymptoticSolve can be used to find these perturbation expansions:

AsymptoticSolve[x^2 - 1 == ϵ x, {x, 1}, {ϵ, 0, 7}]

{{x -> 1 + ϵ/2 + ϵ^2/8 - ϵ^4/128 + ϵ^6/ 1024}}

AsymptoticSolve[x - 2 == ϵ Cosh[x], {x, 2}, {ϵ, 0, 7}]

{{x -> 2 + ϵ Cosh[2] + ϵ^2 Cosh[2] Sinh[2] + 1/2 ϵ^3 (Cosh[2]^3 + 2 Cosh[2] Sinh[2]^2) + 1/3 ϵ^4 (5 Cosh[2]^3 Sinh[2] + 3 Cosh[2] Sinh[2]^3) + 1/24 ϵ^5 (13 Cosh[2]^5 + 88 Cosh[2]^3 Sinh[2]^2 + 24 Cosh[2] Sinh[2]^4) + 1/15 ϵ^6 (47 Cosh[2]^5 Sinh[2] + 100 Cosh[2]^3 Sinh[2]^3 + 15 Cosh[2] Sinh[2]^5) + 1/720 ϵ^7 (541 Cosh[2]^7 + 7746 Cosh[2]^5 Sinh[2]^2 + 7800 Cosh[2]^3 Sinh[2]^4 + 720 Cosh[2] Sinh[2]^6)}}

AsymptoticSolve[x^2 - 1 == ϵ Exp[x], {x, 1}, {ϵ, 0, 7}]
AsymptoticSolve[x^2 - 1 == ϵ Exp[x], {x, -1}, {ϵ, 0, 7}]

{{x -> 1 + (E ϵ)/2 + (E^2 ϵ^2)/8 + (E^3 ϵ^3)/16 + ( 13 E^4 ϵ^4)/384 + (E^5 ϵ^5)/48 + (69 E^6 ϵ^6)/ 5120 + (841 E^7 ϵ^7)/92160}}

{{x -> -1 - ϵ/(2 E) + (3 ϵ^2)/(8 E^2) - (7 ϵ^3)/( 16 E^3) + (235 ϵ^4)/(384 E^4) - (121 ϵ^5)/(128 E^5) + ( 7959 ϵ^6)/(5120 E^6) - (245953 ϵ^7)/(92160 E^7)}}

Answer 4

According to standard perturbation theory for static Hamiltonians of this type, the ground state of the system may be approximated by:

\begin{align} |0\rangle_V &= |0\rangle_H - \epsilon \sum_{\alpha\neq 0} \frac{U_{\alpha 0}}{\alpha} |\alpha\rangle_H \\ \notag &\simeq N^{-\frac12} \left ( |0\rangle_H + \beta_1 |1\rangle_H + \beta_2 |2\rangle_H \right)\,, \end{align} where $N=1+\beta_1^2+\beta_2^2$, $\beta_k = -\epsilon U_{\alpha 0}/\alpha$ and $U_{\alpha 0}={}_H\langle \alpha|U| 0 \rangle_H$, $\alpha=1,2$.

Hope this helps to formulate your problem.