Message Boards Message Boards

0
|
6493 Views
|
22 Replies
|
10 Total Likes
View groups...
Share
Share this post:

trascendental equations in Mathematica

Posted 10 years ago

Hi all,

my name is luca I am a physicist and I do research in a company in Germany. I am trying to get a solution with Mathematica to

Solve[x^a-a*Log[1-x]==c && 1/2<x<1&&0<a<1,x,Reals]

but I cannot get a solution. x is the variable, a and c are parameters. From "plot" I can tell there is one solution. Could anybody suggest the best way to get the root in this case? Solve, NSolve, Reduce and Root would generate the same Warning " This system cannot be solved with the methods available to ..."

Thanks a lot in advance

POSTED BY: Luca D.
22 Replies

my name is luca

I live on the second floor ... ---- sorry, could simply not resist ----.

Have you tried already to do a RegionPlot3D with the equation? This is the surface in question

region plot

This is the companion as surface

contour plot

The goal is a function $x=f(c,a)$ which is as much as possible an analytic expression, right?

P.S.: Originallyl I wrote ParametricPlot3D, which was wrong, of course, I'm sorry about that.

POSTED BY: Udo Krause

You can find numerical solution if you specify values of c and a

In[3]:= c = 1; a = .5; FindRoot[x^a - a*Log[1 - x] == c, {x, .1}]

Out[3]= {x -> 0.468203}
POSTED BY: S M Blinder
POSTED BY: Udo Krause
Posted 10 years ago
POSTED BY: Luca D.
POSTED BY: Henrik Schachner
Posted 10 years ago
POSTED BY: Jim Baldwin
Posted 10 years ago

Hi Henrik,

we'getting there .... thanks for your suggestion :) The Newton's method is working nicely for some values of (x,a). However, sometimes the provided solution is complex and the comparison to the attended values get bad. Would you have suggestion on how to restrict the Newton's method to reals?

POSTED BY: Luca D.

For some purposes this may work:

myRoot[a_?NumericQ, c_?NumericQ] := 
 x /. ToRules@
   Reduce[x^Rationalize[a] - Rationalize[a]*Log[1 - x] == 
      Rationalize[c] && 1/2 < x < 1, x, Reals]

For example, you get a plot with

Plot[myRoot[a, 2], {a, 1/2, 1}]
POSTED BY: Gianluca Gorni

Hi Luca,

I think the complex values show up when x<0. This typically happens when there is no solution at all (e.g. when c>0). Have you ever encountered a complex "solution" when there definitely is a zero crossing?

Ciao Henrik

POSTED BY: Henrik Schachner
Posted 10 years ago

Hi, thanks for your answer. However I rather need an analytical solution (or an approximation of it), as I am going to use the solution at run-time afterwards. Replacing a Taylor expansion of the equation doesn't seem to solve the problem, as Mathematica doesn't output the root: it is clearly one per (x,a) looking at the plot. An idea I had to overcome the problem is to generate sample points x^a-a*Log[1-x] with 1/2<x<1&&0<a<1 and then try to fit it. Replace the function with the fitted function and hopefully I get an equation, which is easier to solve. Well, I don't know if it is the best solution to this problem and actually I cannot do it: I am currently investigation how to get it on the web.

POSTED BY: Luca D.
Posted 10 years ago

Hey I heard already that song :) ... my mistake I shouldn't have said that ;D

Yes I did try with ParametricPlot3D, but it doesn't help. Yup, the goal is x = f(c,a) ...

POSTED BY: Luca D.
Posted 10 years ago

Hi Jim,

thanks for your suggestion. I do need to use the analytical solution outside Mathematica, so I believe "Interpolation" does not help.

POSTED BY: Luca D.
Posted 10 years ago

Hi Gianluca,

thanks for your help. I've just noticed the complex numbers come when I evaluate the Newton's method too far away from the actual solution. This brings the story to the original problem, because I need to know where the solution roughly is before getting an analytical expression of it. Outside Mathematica the solution should be used often, so I might need a different approach to get one analytical expression to be quickly used.

POSTED BY: Luca D.

One searches for a $y$ satisfying the equation $x^a=-a \log(1-y)$

In[16]:= Solve[x^a == -a Log[1 - y] && 0 < a < 1 && 1/2 < x < 1, y, Reals]
Out[16]= {{y -> ConditionalExpression[1 - E^(-(x^a/a)), 0 < a < 1 && 1/2 < x < 1]}}

that's trivial but allows to transform the equation into $(1-x)\exp\left(-\frac{x^a}{a}\right)=\exp\left(-\frac{c}{a}\right)$. This will not be solved again

Solve[(1 - x) Exp[-x^a/a] == Exp[-c/a] && 1/2 < x < 1 && 0 < a < 1, x, Reals]
Solve::nsmet: This system cannot be solved with the methods available to Solve. >>

but can be expressed as $(1-x)\exp\left(\frac{c-x^a}{a}\right)=1$.

It's interesting to see how different quality the corresponding ContourPlot3D for the three equivalent equations are.

Any derivative of the above mentioned expression must vanish, but I do not see an usage of that fact.

But one could try to develop around $a = 1$. At least for $a=1$ a solution appears

In[36]:= Solve[(1 - x) Exp[c - x] == 1 , x]
During evaluation of In[36]:= Solve::ifun: Inverse functions are being used by Solve, ome solutions may not be found; use Reduce for complete solution information. >>
Out[36]= {{x -> 1 - ProductLog[E^(1 - c)]}}
POSTED BY: Udo Krause
POSTED BY: Udo Krause
Posted 10 years ago
POSTED BY: Jim Baldwin
Posted 10 years ago

I found an approach at https://www.physicsforums.com/threads/how-to-solve-log-x-x-2-0-for-x.537411/ which uses a fixed point iterative solution.

For more numeric stability I rearranged x^a - a Log[1-x] == c to x == 1 - Exp[(x^a - c)/a]. Using a starting value of 1:

f[a_, c_, n_: 20] := Module[{z, i},
  (* n is number of iterations *)
  z = 1;
  Do[
   z = 1 - Exp[(z^a - c)/a],
   {i, n}];
  z]

(* Example *)
a = 0.5;
c = 10;
x = f[a, c]

(* Check:  the value below should be close to zero *)
If[x < 1, x^a - a Log[1 - x] - c, "Estimate is too close to 1"]

(* Use FindRoot now that we hopefully have a good starting value *)
FindRoot[z^a - a Log[1 - z] == c, {z, x}]

When the desired result is too close to 1, then things don't work so well. If you get a solution less than 0.9999, you're probably safe. Using variables with higher than double precision would also help. Again, because of the linear relationship of a and c for a specific value of x, one could determine if the result will be bigger than some safe threshold like maybe 0.9999.

If you need to know if x is something like 0.99999999994, then you'd absolutely need multiple precision arithmetic for your "outside" program.

POSTED BY: Jim Baldwin

This returns complex $x$

f[a, c, n_: 20] := Module[{z, i},(*n is number of iterations*)z = 1; Do[z = 1 - Exp[(z^a - c)/a], {i, n}]; z]

as it stands

In[21]:= x = f[0.5, 0.7, 10000]
Out[21]= 1.29175 - 0.278866 I

In[22]:= x^0.5 - 0.5 Log[1 - x] - 0.7
Out[22]= 0.896757 - 1.31137 I

not belonging to 1/2 < x < 1.

POSTED BY: Udo Krause

Let's set a=1/2, because of $1=(1-x)\exp\left(\frac{c-x^a}{a}\right)$ one could ask for the fourier transform. The 1 goes to the dirac delta , but

In[64]:= (* a = 1/2 *)
FourierTransform[E^(2 (c - Sqrt[x])) (1 - x), x, w]

Out[64]= (1/(8 Sqrt[2 \[Pi]] w^4 Abs[w]^(9/2)))E^(
 2 c - I/w) (-2 I E^(2 I/w) Sqrt[2 \[Pi]] w^5 (2 - 3 I w + 2 w^2) + 
   Sqrt[2 \[Pi]]
     w^4 (-3 I w + E^(2 I/w) (4 - 3 I w + 4 w^2)) Abs[w] - 
   3 I (-1 + E^(2 I/w)) Sqrt[2 \[Pi]] w^3 Abs[w]^3 + 
   4 E^(I/w) Abs[w]^(
    9/2) (E^(I/w) Sqrt[\[Pi]] Sqrt[I w] (2 I + 3 w + 2 I w^2) - 
      2 I w (1 - I w + w^2) + 
      E^(I/w) Sqrt[\[Pi]] Sqrt[
       I w] (2 - 3 I w + 2 w^2) Erfi[1/Sqrt[I w]]) + 
   32 I E^(I/w) w^5 Sqrt[Abs[w]]
     HypergeometricPFQ[{2}, {3/4, 5/4}, -(1/(4 w^2))] + 
   8 E^(I/w) w^4 Abs[w]^(
    5/2) (I w HypergeometricPFQ[{1}, {1/4, 3/4}, -(1/(4 w^2))] - 
      2 HypergeometricPFQ[{1}, {3/4, 5/4}, -(1/(4 w^2))] + 
      HypergeometricPFQ[{1, 3/2}, {1/4, 1/2, 3/4}, -(1/(4 w^2))]))

no chance to solve for w and to transform back. You need a truly clever trick to transform the condition or a numerical procedure along the lines Jim Baldwin proposed.

POSTED BY: Udo Krause
Posted 10 years ago

Ooops! Not sure where I got that c must be greater than 1.19315 as c must be greater than or equal to 1. However, the numerical method described above does give a real result when a is between 0 and 1 and c is greater than zero.

POSTED BY: Jim Baldwin
Posted 10 years ago
POSTED BY: Jim Baldwin
Posted 10 years ago

Sorry for one more post. c must be greater than 0.5^a - a*Log[0.5] so that 0.5 < x < 1. And rather than fixing a number of iterations, it is almost certainly better to have a minimum and maximum number of iterations along with a reasonable stopping rule. I was just being lazy with the example.

POSTED BY: Jim Baldwin
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract