How do I find a series solution for $e^{-\frac{1}{2}f'(x)} \mathrm{cosh}( f(x) ) = ax + b$?

Question

I am trying to approximate a function $f(x)$ satisfying a relation between $f(x)$ and its first derivative. How do I find a series solution for $$e^{-\frac{1}{2}f'(x)} \mathrm{cosh}( f(x) ) = ax + e^{-\frac{1}{2}f'(0)}$$ about $x=0$? Thus far I've looked at the numerical solution.

s = NDSolve[ { -0.5*y'[x] + Log[Cosh[y[x]]] == 
    1 + Log[ 0.1*y[x] + Exp[-1 - 0.25 ]] , y[0] == 0}, y, {x, 0, 4} ]
Plot[Evaluate[ y[x] /. s], {x, 0, 4}, PlotRange -> All ]

And if this were a linear ODE I could use Series[DifferentialRoot[ ] ].

Answer 1

Series Expansion

Answering the question we show how to construct the series expansion of f[x] about x=0.

From the basic equation we find

f'[x]==2(Log[a x + b]-Log[Cosh[f[x]]])

the derivative of f[x] in terms of the function f[x] itself and some other terms.

Hence in the Taylor expansion of the function

f[x] = f[0] + x f'[0] + x^2/2! f''[0] + x^3/3! f'''[x] + ...

we can replace f'[x] by

g[x_] = 2 (Log[a x + b] - Log[Cosh[f[x]]]);

f''[x] is then replaced by g'[x] which in turn is expressed by f[x] via

g'[x]

(* Out[62]= a/(b + a x) - Tanh[f[x]] Derivative[1][f][x] *)

where f'[x] on the right hand side must be replaced by g[x] and so on.

This means that all coefficients of the Tayor expansion are expressed in terms of f[0], a, and b.

In order to formalize this procedure in Mathematica we define recursively the auxiliary functions

g[k_, x_] := D[g[k - 1, x], x] /. f'[x] -> g[x]
g[1, x_] = g[x];

The first few functions are

Table[g[k, x], {k, 1, 3}] // Column

$\begin{array}{l} \{\text{Log}[b+a x]-\text{Log}[\text{Cosh}[f[x]]]\} \\ \left\{\frac{a}{b+a x}-(\text{Log}[b+a x]-\text{Log}[\text{Cosh}[f[x]]]) \text{Tanh}[f[x]]\right\} \\ \left\{-\frac{a^2}{(b+a x)^2}-(\text{Log}[b+a x]-\text{Log}[\text{Cosh}[f[x]]])^2 \text{Sech}[f[x]]^2-\text{Tanh}[f[x]] \left(\frac{a}{b+a x}-(\text{Log}[b+a x]-\text{Log}[\text{Cosh}[f[x]]]) \text{Tanh}[f[x]]\right)\right\} \\ \end{array}$

Hence the requested series expansion to order n is given by

fs[x_, n_] := f[0] + Sum[x^k/k! (g[k, x] /. x -> 0), {k, 1, n}]

Example n=3

fs[x, 3]

$\text{f0}+x (\text{Log}[b]-\text{Log}[\text{Cosh}[\text{f0}]])+\frac{1}{2} x^2 \left(\frac{a}{b}-(\text{Log}[b]-\text{Log}[\text{Cosh}[\text{f0}]]) \text{Tanh}[\text{f0}]\right)+\frac{1}{6} x^3 \left(-\frac{a^2}{b^2}-(\text{Log}[b]-\text{Log}[\text{Cosh}[\text{f0}]])^2 \text{Sech}[\text{f0}]^2-\text{Tanh}[\text{f0}] \left(\frac{a}{b}-(\text{Log}[b]-\text{Log}[\text{Cosh}[\text{f0}]]) \text{Tanh}[\text{f0}]\right)\right)$

Numerical calculation of f(x) via an ODE

This was my original text, which, however, is not an answer to the problem.

Solving for f' you can write your equation as an ODE

eq = f'[x]/2 == Log[a x + b] - Log[Cosh[f[x]]];

Then solve it numerically after having defined the values of the parameters a and b.

For example

a = 1; b = 1;
ff[x_] = f[x] /. 
  NDSolve[f'[x]/2 == Log[a x + b] - Log[Cosh[f[x]]] && f[0] == 1, f[x], {x, -2, 10}][[1]]

Plot the function

Plot[ff[x], {x, -2, 10}, PlotRange -> {0, 10}]
(* 141107_NDSolve.jpg *)

You can play with the values of the parameters and the initial condition to adapt the problem to the situation you want.

Regards,
Wolfgang

Answer 2

I stumbled across this question just now, but it seems to me that there's a much easier way to do this problem via user-defined SeriesData objects. If we want to define $$ f(x) = c_0 + c_1 x + c_2 x^2 + \dots $$ it's not too hard to do this in Mathematica:

n = 4;  
f[x_] = SeriesData[x, 0, Array[c, n, 0], 0, n, 1]

(* c[0] + c[1] x + c[2] x^2 + c[3] x^3 + O[x]^4 *)

Once we've done this, it's pretty easy to get a set of rules relating the $c_i$ coefficients and $a$:

SolveAlways[Exp[-f'[x]/2] Cosh[f[x]] == a x + Exp[-f'[0]/2], x];     
FullSimplify[%, Assumptions -> Element[c[_], Reals]]

(* {{a -> E^(-(c[1]/2)) (-c[2] + c[1] Sinh[c[0]]),
     c[3] -> 1/3 (c[1]^2 + c[2]^2 - 2 (-1 + c[1]) c[2] Sinh[c[0]]), 
     Cosh[c[0]] -> 1}} *)

Disentangling the above equations, this implies that $c_0 = 0$, $c_2 = -a e^{c_1/2}$, and $c_3 = \frac{1}{3} (c_1^2 + a^2 e^{c_1})$. The coefficient $c_1$ is freely specifiable. Note that we must have $c_0 = 0$: when we plug $x = 0$ into the original equation, we get $\cosh(f(0)) = 1$, or $f(0) = 0$.

If you want to get a direct solution for the $c_i$ coefficients, that's possible too; it just takes a little more work:

Thread[CoefficientList[Normal[Exp[-f'[x]/2] Cosh[f[x]] - (a x + Exp[-f'[0]/2])], x] == 0];
Solve[%, Drop[Array[c, n, 0], {2}]]

(* {{c[0] -> 0, c[2] -> -a E^(c[1]/2), c[3] -> 1/3 (a^2 E^c[1] + c[1]^2)}} *)

The Drop is necessary so Mathematica doesn't try to solve for the free coefficient $c_1$. (There's probably a nicer way to do that.) This method has the advantage that it keeps working even at higher values of $n$; my first method conks out around $n = 9$ (complaining about the relations being "essentially non-algebraic"), but the second method spits out an answer in about a second.

Answer 3

Using new function AsymptoticDSolveValue form Mathematica 11.3:

sol = AsymptoticDSolveValue[{-y'[x]/2 + Log[Cosh[y[x]]] == 
1 + Log[y[x]/10 + Exp[-(5/4)]], y[0] == 0}, y[x], {x, 0, 6}]

(* x/2 - 1/20 E^(5/4) x^2 + 
1/240 (20 + E^(5/2)) x^3 + ((-800 E^(5/4) - 
17 E^(15/4)) x^4)/48000 + ((70000 + 10400 E^(5/2) + 
159 E^5) x^5)/4800000 + ((-570000 E^(5/4) - 35900 E^(15/4) - 
487 E^(25/4)) x^6)/144000000*)

Comparison with numeric solution:

numeric = NDSolve[{-0.5*y'[x] + Log[Cosh[y[x]]] == 
1 + Log[0.1*y[x] + Exp[-1 - 0.25]], y[0] == 0}, y, {x, 0, 1}]

Plot[{sol, Evaluate[y[x] /. numeric]}, {x, 0, 1}, 
PlotLegends -> {"Series", "Numeric"}]