How to solve ArcSinh equation for many differents parameters?

I am totally stuck trying to solve the equation below

(1/2) x Sqrt[1 + 4 a^2 x^2] + ArcSinh[2 a x]/(4 a) == d

Where a and d are given numbers (a is always equal to 0.022). The thing is that I would like to obtain x for different values of parameter d, stored in a .TXT file

Please, do anyone know how to tell Mathematica to give me the different x's for each d from the file? How to export the resulting vector to excel? I have already tried with NSolve and Solve, but nothing happens.

• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Read the faq! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign! Oct 30 '14 at 18:44

Introduction

There are a few ways to make such a table: (1) a symbolic way that seems conceptually clear and whose slowness is not prohibitive for a table of at most a few thousand entries; (2) using the table created by NDSolve in solving the equation by integrating its derivative either (a) directly or (b) correcting the errors in it; and (3) making a regular-interval table in two ways (my feeling is this should be unnecessary and it's a hold-over from human-use, printed tables). See also my answer to https://mathematica.stackexchange.com/a/27984, where I did some of these things.

Basically, we want a table of the form

Table[{d, x}, {d,...}]

where d = f[x] is given and x = InverseFunction[f][d] is sought. One of the simplest approaches is to reverse the entries in a table of the form

Table[{x, d}, {x,...}]

Finally, an equation f[x[d]] == d may differentiated to get f'[x[d]] x'[d] == 1 and solved with NDSolve, given that we know or can find an initial value. These relationships between the variables x and d may be applied to any equation f[x, d] == 0, except InverseFunction[f][g[d]] would apply only to equations in which the variables can be separated, f[x] == g[d], and InverseFunction[f] is computable (c.f. the inverse function theorem).

A table has to have a range. Since none was indicated, I'll pick 0 <= d <= maxd, where maxd = 100. The OP's function has odd symmetry, so it is unnecessary to create a table for negative values of d.

Some definitions

I'll keep the definitions in terms of a symbolic a. One can use Block[{a}, code] to keep it symbolic through part of a computation. It is occasionally helpful.

Needs["GeneralUtilities`"];

timing = AccurateTiming;  (* Use AbsoluteTiming in V9 or search this site for timeAvg *)

a = 0.022;
lhs := (1/2) x Sqrt[1 + 4 a^2 x^2] + ArcSinh[2 a x]/(4 a)

dfn = Block[{a}, Function @@ {{x}, lhs}]  (* d as a function of x *)
(*  Function[{x}, 1/2 x Sqrt[1 + 4 a^2 x^2] + ArcSinh[2 a x]/(4 a)]  *)

Block[{a}, dprime := Evaluate@Simplify[dfn'[x], x >= 0]];
Definition[dprime] maxd = 100;
maxx = x /. FindRoot[lhs == maxd, {x, 0}]
(*  58.5729  *)

Initial condition: d equals zero at x equals zero.

dfn
(*  0.  *)

Symbolic solution

DSolve can be used to construct the InverseFuction that defines x as a function of d. (DSolve works better with a symbolic a than with an approximate 0.022.)

xfnD = Block[{a}, DSolveValue[{D[dfn[x[d]], d] == 1, x == 0}, x, d]]
(*  Function[{d},
InverseFunction[ArcSinh[2 a #1]/(4 a) + 1/2 #1 Sqrt[1 + 4 a^2 #1^2] &][d]]  *)

NDSolve solution

One small trick: We need to change x to a function x[d] in the derivative dprime. We're using the inverse function theorem, $x'(d) = 1/d'(x)$, if I may abuse notation, to construct x as an interpolating function of d. NDSolve has many options to control how the integration is done, including the step size and the precision of the result..

xfnN = NDSolveValue[{x'[d] == 1/(dprime /. x -> x[d]), x == 0},
x, {d, 0, maxd}, MaxStepSize -> 0.1];

I want to compare the timings of the methods, so I'll save it:

ndsolvetime = timing[
NDSolveValue[{x'[d] == 1/(dprime /. x -> x[d]), x == 0}, x, {d, 0, maxd},
MaxStepSize -> 0.1];
]
(*  0.00448259  *)

Making tables

We can make tables of values of {d, x} either in the form

tab = Table[{d, xfn[d]}, {d,...}]

or

tab = Table[{lhs, x}, {x,...}]

If we have a table tab, we can create an (order-3) interpolating function with

xTest = Interpolation[tab]

or, if the tables are to be used for linear interpolation,

xTest = Interpolation[tab, InterpolationOrder -> 1]

One might be concerned with the accuracy of the table. We can estimate the maximum relative error of the approximation with the following function. Sometimes it gives a warning about convergence failure, it's hardly important for this purpose. A new seed is used each time, so one can reevaluate to examine the stability of the estimate. The relative error tends to grow near x == 0 although the absolute error diminishes. Consequently, we'll maximize the error over 1 < x < maxx by default. (Note that using cubic interpolation generally improves the error quite significantly, but I'll test the linear interpolation in what follows.)

maxrelerr[xfn_, xmin_: 1] :=
NMaximize[{Abs[(xfn[dfn[x]] - x)/x], xmin < x < maxx}, x,
Method -> {"RandomSearch", "SearchPoints" -> 300, "RandomSeed" -> RandomInteger}]

Easiest: regular intervals of x

This method has the advantage of being exact (within machine precision, that is) on the points in the table.

nentries = 1000;
tab = Table[{lhs, x}, {x, 0., maxx, maxx/(nentries - 1)}]; // timing
(*  0.00625478  *)

xTest = Interpolation[tab, InterpolationOrder -> 1];
maxrelerr[xTest]
(*  {8.30295*10^-7, {x -> 1.02577}}  *)

Using the inverse function: regular intervals of d

This is also easy, although it is a bit slow.

tab = Table[{d, xfnD[d]}, {d, 0, maxd, 1/10}]; // timing
(*  0.710307  *)

xTest = Interpolation[tab, InterpolationOrder -> 1];
maxrelerr[xTest]
(*  {2.41031*10^-6, {x -> 1.04883}}  *)

Interestingly the first method seems more accurate. One should keep in mind that it will depend on the function being interpolated. This one is nice and gets straighter as d increases. The previous method happens to take larger steps, for this particular function, as d increases. It won't always happen that way.

This method is similar to Mr.Wizard's, albeit this one is a little slower.

Table[Quiet@FindRoot[lhs == d, {x, 0}], {d, 0, maxd, 1/10}]; // timing
(*  0.479001  *)

Extracting the NDSolve table

The table is stored inside the InterpolatingFunction (see What's inside InterpolatingFunction[{{1., 4.}}, <>]? for more). The list of values of d is given by

First @ xfnN["Coordinates"]

and the list of values of x is given by

xfnN["ValuesOnGrid"]

We can form the table via transposition of the lists.

tab = Transpose[{First@xfnN["Coordinates"], xfnN["ValuesOnGrid"]}]; // timing
tab // Length
(*  0.000373258
1008         *)

The total timing is roughly

ndsolvetime + 0.000373258
(*  0.00485584  *)

xTest = Interpolation[tab, InterpolationOrder -> 1];
maxrelerr[xTest]
(*  {2.40168*10^-6, {x -> 1.04722}}  *)

There is not really an advantage to using NDSolve in this case for two reasons. First, the simple form of the equation f[x] == d allows the first method to work. Second, the lhs is monotonic and the second derivative monotonically tends to zero. This method is more useful for more complicated equations.

Aside from interpolation, one source of error in this method is that NDSolve calculates the table entries only approximately (with a precision of around one half of machine precision). These usually can be improved to machine precision by a single step of Newton's method, which I'll share:

newx = x - (lhs - d)/dprime /.
{d -> First@xfnN["Coordinates"], x -> xfnN["ValuesOnGrid"]};
newtab = Transpose[{First@xfnN["Coordinates"], newx}];

xTest1 = Interpolation[tab];       (* cubic interpolation *)
xTest2 = Interpolation[newtab];

{maxrelerr[xTest1] // Quiet,
maxrelerr[xTest2] // Quiet}

{maxrelerr[xTest1, 20] // Quiet,   (* error for larger d *)
maxrelerr[xTest2, 20] // Quiet}

(*  {{9.07249*10^-9,  {x -> 1.00002}}, {1.12738*10^-10, {x -> 1.04673}}}
{{3.41066*10^-10, {x -> 20.0006}}, {8.87196*10^-13, {x -> 26.5654}}}  *)

Again since the graph is straightening out as d gets larger, the relative error is improved greatly for larger values of d because the error is primarily due to the error in the table.

Using NDSolve to construct regular intervals in d

tab = Table[{d, xfn[d]}, {d, 0, 100, 1/10}]; // timing
ndsolvetime + %
(*  0.00265543
0.00713802 *)

xTest = Interpolation[tab, InterpolationOrder -> 1];
maxrelerr[xTest]
(*  {2.38447*10^-6, {x -> 1.04884}}  *)

The relative error is not quite as good as the first method but comparable to the others.

Exporting

To export to Excel is easy:

Export[FileNameJoin[{\$TemporaryDirectory, "foo.xls"}], tab]
(*  "/var/folders/9d/68khy4s15sjf9qfpnhqz9tnc0000gn/T/foo.xls"  *)

SystemOpen[%]

Other formats may be found in the Listing of All Formats.

Summary

All the methods are reasonably fast for a one-time table generation. The accuracy of each may be controlled by controlling the step size. I think it is important to examine the accuracy of the generated table and adjust the steps to get the desired accuracy. For the OP's particular example, the first method wins hands-down. NDSolve is not shown to advantage in this problem, but it is normally what I use for getting a numerical approximation to a function.

• Bored, Michael? :-O Oct 30 '14 at 21:35
• @Mr.Wizard Embarrassingly, no. This has come up before, at least the way I look at it. I've solved it several times in several ways myself, sometimes for myself. (I have a way of parametrizing an implicit curve by arc length, but I have yet to find a question to use it in an answer. This is close, but they wanted y = f[x].) Here I tried to give an answer to a generalized question and hope future site visitors might find it useful. I asked two of the chief numerics guys are WRI about this kind of problem, and so I'm feeling confident, too. Oct 30 '14 at 22:53

I'll try to address this one step at a time. Since a is assigned a value of 0.022 I shall assume we will be working numerically. Let's plot our function.

a = 0.022;

expr = (1/2) x Sqrt[1 + 4 a^2 x^2] + ArcSinh[2 a x]/(4 a)

Plot[expr, {x, -1*^5, 1*^5}]
1/2 x Sqrt[1 + 0.001936 x^2] + 11.3636 ArcSinh[0.044 x] There appear to be no complications so this should be straightforward. We can use FindRoot to search for a solution:

FindRoot[expr == 1318.77, {x, 0}]
{x -> 241.055}

You can find a series of solutions using Function and Map

values = {37, 192.5, 281.333, 9488.53};

sols = FindRoot[expr == #, {x, 0}] & /@ values
{{x -> 29.8589}, {x -> 86.2576}, {x -> 106.631}, {x -> 654.939}}

Or using Table:

sols = Table[FindRoot[expr == d, {x, 0}], {d, values}]
{{x -> 29.8589}, {x -> 86.2576}, {x -> 106.631}, {x -> 654.939}}

You can convert this result to simple numbers using ReplaceAll:

x /. sols
{29.8589, 86.2576, 106.631, 654.939}

There are already lots of examples of importing from and exporting to files on this site, as well as the documentation for Import and Export so I will not address that here unless you are experiencing specific difficulties.