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

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[ ] ].

• Could you please post the Mathematica code for your equations so we don't have to retype them? – Dr. belisarius Nov 7 '14 at 5:24
• I'm entirely new to Mathematica. I looked up how to define differential operators and nobody seems to agree. So I'm looking for best practice. – Quant Nov 7 '14 at 5:26
• So I recommend searching this site for DSolve[] before asking – Dr. belisarius Nov 7 '14 at 5:41
• Related: 25363, 26938 – Michael E2 Nov 7 '14 at 11:36
• (1) Always helpful to write the input as explicit Mathematica (InputForm) code. Something like Exp[-1/2*D[f[x],x]]*.... (2) Possibly what you want is to subtract rhs from lhs, take ser = Series[..., {x,0,4}] for example, and then SolveAlways[ser,x]. (I'd test it but you didn't include that input...) – Daniel Lichtblau Nov 7 '14 at 15:46

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

• Another way to accomplish the effect of g: dfrule = f'[x] -> 2 (Log[a x + b] - Log[Cosh[f[x]]]); NestList[D[First[#], x] -> (D[Last[#], x] /. dfrule) &, dfrule, 3] /. x -> 0, where 3 may be replaced by any k. – Michael E2 Nov 7 '14 at 22:12
• You've not only answered my question but got me up to speed with the software. Thanks so much. – Quant Nov 7 '14 at 22:17
• Thanks to the (unknown) upvoters. – Dr. Wolfgang Hintze Jun 26 '15 at 11:53

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.

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"}]