# trascendental equations in Mathematica

Posted 8 years ago
5590 Views
|
22 Replies
|
10 Total Likes
|
 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
22 Replies
Sort By:
Posted 8 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 8 years ago
 Now I remember. When a = 1, then c must be greater than 1.19315 to have 0.5 < x < 1. So while the numeric method provides real solutions for c > 1, some of those will result in x <= 0.5. A test could be added to check for the proper lower limit of c given a.
Posted 8 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 8 years ago
 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 8 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 8 years ago
 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 8 years ago
 I think the solution for you (as you need a good algorithm for the outside world - C? Fortran? Basic? Perl?) will ultimately depend on what values of c and a you expect and how precise you want to be able to determine x.From everything folks have mentioned above: c must be greater than or equal to 1.19315 (as 0.5 < x < 1). For any value of c the minimum value of x will be when a = 1. This will give you a lower limit for x in any search algorithm. You can use what Udo Krause mentioned to find that lower limit. If you only need to know x up to 0.9999, then you won't need to perform any iterative method if c is greater than about 10. You just return 0.9999 or 1.0000. Also, because c is linearly related to a for a given value of x, if there is some maximum x like 0.9999, then one can determine if the (c,a) pair is above the line for x = 0.9999 which the function could then return 0.9999 if that was your upper threshold for x. Depending on the precision desired for x, you might need to consider multiple precision support as double precision might not give enough precision. When c = 100, then the maximum value of x is 0.99999999999999999999999999999999999999999989887785073895514700547142\ 0237974134503800561929733764673745574785011843276329550445380987673228\ 2904164154378408341327050646705054851432685561291692934129708640969732\ 9364105768971523513419413703741965175554141229402839267590651808760914\ 99111583603367054544515.. I don't think that double-precision will get you there (unless maybe there's merit in finding 1 - x).
Posted 8 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 8 years ago
 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 8 years ago
 I know you want an analytic expression but if all of your future calculations are within Mathematica, why not use interpolation functions? (Or, is the need for an analytic expression so that you can use the results outside of Mathematica?) n = 100; (* Number of elements in each row of the interpolation table \ *) (* Set ranges for a and x *) amin = 0.; amax = 1.; xmin = 0.5; xmax = 1.; (* Create table of values of a, c, and x *) t = Flatten[ Table[{{amin + (amax - amin) i/ n, (xmin + (xmax - xmin) j/n)^(amin + (amax - amin) i/ n) - (amin + (amax - amin) i/n)* Log[1 - (xmin + (xmax - xmin) j/n)]}, xmin + (xmax - xmin) j/n}, {i, 1, n}, {j, 1, n - 1}], 1]; (* Create interpolation function *) f = Interpolation[t, InterpolationOrder -> 1] (* Contour plot *) ContourPlot[f[a, c], {a, 0.01, 0.99}, {c, 1.01, 6}] (* Any particular point: f[a,c] *) f[0.75, 2]  Out[672]: 0.788073 
Posted 8 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 8 years ago
 I tried Simplify[x^a-a*Log[1-x]-c==0,1/2<=x<=1&&Element[x,Reals]&&0
Posted 8 years ago
 Hey I heard already that song :) ... my mistake I shouldn't have said that ;DYes I did try with ParametricPlot3D, but it doesn't help. Yup, the goal is x = f(c,a) ...
Posted 8 years ago
 Can't you represent the power $x^a$ as a logarithmic expression to homogenize the look of the equation? Have you tried all flavours of simplify on the equation, to get some idea where to find the next step?
Posted 8 years ago
 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, so some solutions may not be found; use Reduce for complete solution information. >> Out[36]= {{x -> 1 - ProductLog[E^(1 - c)]}} 
Posted 8 years ago
 For $a \rightarrow 0$ you catch a condition: $c \rightarrow 1$. Somehow raising $a$ from $0$ one could possibly see what $c$ has to be? In[44]:= Solve[E^((c - x^a)/a) (1 - x) == 1, c] Out[44]= {{c -> ConditionalExpression[x^a - a (2 I \[Pi] C[1] + Log[1 - x]), C[1] \[Element] Integers]}} Returned to sender. Gosh! Solve[E^((c - x^a)/a) (1 - x) == 1, a] Solve::nsmet: This system cannot be solved with the methods available to Solve. >> of course.You could solve In[47]:= Solve[E^((c - x^a)/a) (c - x^a)/a == 1, x] During evaluation of In[47]:= Solve::ifun: Inverse functions are being used by Solve, so some solutions may not be found; use Reduce for complete solution information. >> Out[47]= {{x -> (c - a ProductLog[1])^(1/a)}} but again no benefit from this information is visible.
Posted 8 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
Posted 8 years ago
 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 questionThis is the companion as surfaceThe 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 8 years ago
 Hi Luca,maybe this simple (an hopefully not too stupid) idea is helpful: Using Newtons method with analytic expressions!Define your function and the Newton iteration step: f[x_] := x^a Log[1 - x] - c newton[x_] := x - f[x]/f'[x] Then you do some (not too many, e.g. 3) interations: sol[x_] = Nest[newton, x, 3]; This results in a large/huge - but analytical ! - formula. Since you decided to have 1/2 < x < 1 you can start with an initial value for e.g. x=0.75; an approximation of the solution then should be (depending on the parameters a and c) as an example: sol[.75] /. {a -> .5, c -> -.7} Greetings Henrik
Posted 8 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 8 years ago
 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 8 years ago
 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} `
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments