# Power series solution for ODE

The ODE I have is $$y'(x)+e^{y(x)}+\frac{e^x-e^{-x}}{4}=0, \hspace{0.2cm} y(0)=0$$ I want to determine the first five terms (coefficients $a_0,\ldots, a_5$) of the power series solution $$y(x)=\sum_{k=0}^{\infty} a_kx^k$$ So far, I know that $$y'(x)=\sum_{k=1}^{\infty} a_kkx^{k-1}$$ Now I plug these back into the equation and get: $$\sum_{k=1}^{\infty} a_kkx^{k-1} + e^{\sum_{k=0}^{\infty} a_kx^k} + \frac{e^x-e^{-x}}{4}=0$$. Now I'm not sure how to continue with this. Please help.

• I didn't know Lebron was studying ODE, but you might have a better time by expanding $e^y$ in terms of its power series, then plug in the power series approximation for $y$ from there. – DaveNine Mar 27 '16 at 22:08
• Could you elaborate on that a bit? – LeBron James Mar 28 '16 at 2:45
• $e^y = \sum_{k=0}^{\infty} \frac{y^k}{k!}$, use this to make computing the coefficients easier. I'd consider truncating the series in order to compute the coefficients you want as well. – DaveNine Mar 28 '16 at 3:35
• I used the expansion for $e^x$ and $e^{-x}$, but I fail to see how this will make $e^y$ easier, now I have $$e^y = \sum_{k=0}^{\infty}\frac{\sum_{k=0}^{\infty} a_kx^k}{k!}$$ – LeBron James Mar 28 '16 at 3:55
• Another hint, you're going to want to move over the taylor series of $\frac{1}{2} \cosh(x)$ term to the other side, and then match coefficients. This should give you a means to solve for coefficients. – DaveNine Mar 28 '16 at 3:56

## 4 Answers

Since you just need a few terms, setting $$y=\sum_{n=1}^6 a_i x^i$$ (because of the condition $y(0)=0$) you could develop $e^y$ as a Taylor series around $x=0$ and get $$e^y=1+a_1 x+\frac{1}{2} \left(a_1^2+2 a_2\right) x^2+\frac{1}{6} \left(a_1^3+6 a_2 a_1+6 a_3\right) x^3+$$ $$\frac{1}{24} \left(a_1^4+12 a_2 a_1^2+24 a_3 a_1+12 a_2^2+24 a_4\right) x^4+$$ $$\frac{1}{120} \left(a_1^5+20 a_2 a_1^3+60 a_3 a_1^2+60 a_2^2 a_1+120 a_4 a_1+120 a_2 a_3+120 a_5\right) x^5+$$ $$\frac{1}{720} \left(a_1^6+30 a_2 a_1^4+120 a_3 a_1^3+180 a_2^2 a_1^2+360 a_4 a_1^2+720 a_2 a_3 a_1+720 a_5 a_1+120 a_2^3+360 a_3^2+720 a_2 a_4+720 a_6\right) x^6$$

Expanding the term $\frac{e^x-e^{-x}}{4}$ as a Taylor series too, the differential equation then write $$(a_1+1)+\left(a_1+2 a_2+\frac{1}{2}\right) x+\frac{1}{2} \left(a_1^2+2 a_2+6 a_3\right) x^2+\frac{1}{12} \left(2 a_1^3+12 a_2 a_1+12 a_3+48 a_4+1\right) x^3+$$ $$\frac{1}{24} \left(a_1^4+12 a_2 a_1^2+24 a_3 a_1+12 a_2^2+24 a_4+120 a_5\right) x^4+$$ $$\frac{1}{240} \left(2 a_1^5+40 a_2 a_1^3+120 a_3 a_1^2+120 a_2^2 a_1+240 a_4 a_1+240 a_2 a_3+240 a_5+1440 a_6+1\right) x^5+$$ $$\frac{1}{720} \left(a_1^6+30 a_2 a_1^4+120 a_3 a_1^3+180 a_2^2 a_1^2+360 a_4 a_1^2+720 a_2 a_3 a_1+720 a_5 a_1+120 a_2^3+360 a_3^2+720 a_2 a_4+720 a_6\right) x^6=0$$

Cancelling the coefficients lead to $$a_1=-1\qquad a_2=\frac{1}{4}\qquad a_3=-\frac{1}{4}\qquad a_4=\frac{7}{48}\qquad a_5=-\frac{19}{160}$$

I hope I did not make any mistake since my results do not coincide with DaveNine's answer.

• The issue with this is that, through $e^y$, each coefficient $a_1,a_2,...$ are expressed as an infinite series themselves. This means your answers are more of an approximation rather than exact. – DaveNine Mar 28 '16 at 5:50
• @DaveNine. I totally agree with you ! I must confess that I was not very comfortable after writing. Your approach is definitely more rigorous. What I do not understand is why, even for the very first terms, I do not get the same coefficients as yours. Cheers. – Claude Leibovici Mar 28 '16 at 5:55
• Ah I see. Well each coefficient can be represented as a power series because of the exponential term, that's probably why. If you truncate that series then you approximate every coefficient no matter the location in the series of $y$. – DaveNine Mar 28 '16 at 6:07
• thanks, I got the same answers but using a repetitive derivative method, basically solving for $y'(0), y'(0), y''(0)...etc$, then using the Taylor formula to get the coefficients. – LeBron James Mar 28 '16 at 6:17
• I was incorrect, I expanded the non-derivative term incorrectly (You don't even need to anyways.) – DaveNine Mar 28 '16 at 8:05

The answers provided are excellent, but I'll offer what I think is the easiest solution: Just take derivatives of the equation, plug in $0$ and get as many $y^{(n)}(0)$'s as you want, then finally construct your taylor series:

$$y=\sum_{n=0}^\infty \frac{y^{(n)}(0)}{n!}x^n.$$

$$y'+e^{y}+\frac{1}{2}\sinh(x)=0, \quad\quad y(0)=0$$ Plugging in $y(0)=0$, we get $y'(0)+1=0$, thus $\boxed{y'(0)=-1}$.

Now differentiate the original equation as many times as needed to get the number of coefficients that you want: $$y''+y'e^y+\frac{1}{2}\cosh(x)=0$$ $$y'''+(y')^2e^y+y''e^y+\frac{1}{2}\sinh(x)=0$$ etc...

Plugging in $x=0$ into the above equations gives: $$y''(0)+y'(0)e^{y(0)}+\frac{1}{2}\cosh(0)=0 \quad \Rightarrow \quad \boxed{y''(0)=\frac{1}{2}}$$ $$y'''(0)+(y'(0))^2e^y+y''(0)e^{y(0)}+\frac{1}{2}\sinh(0)=0 \quad \Rightarrow \quad \boxed{y'''(0)=-\frac{3}{2}}$$

Thus: \begin{aligned} y&=y(0)+y'(0)x+\frac{1}{2!}y''(0)x^2+\frac{1}{3!}y'''(0)x^3+\cdots \\ &=-x+\frac{1}{2}\frac{1}{2!}x^2-\frac{3}{2}\frac{1}{3!}x^3+\cdots \\ &=-x+\frac{1}{4}x^2-\frac{1}{4}x^3+\cdots \end{aligned}

This method is nice because you don't have to work out recursion relations, etc.

To elaborate on my hint (it my not be the best method, but it should work.)

We have

$$\sum_{k=1}^{\infty} ka_{k}x^{k-1} + \exp \left( \sum_{n=0}^{\infty} a_n x^n \right) = \frac{1}{4}(\exp{(x)} - \exp{(-x)})$$

Setting $x = 0$ and using $a_0 = 0$ we've

$$1+a_1 = 0$$

which gives $a_1 = -1$.

To get other coefficients, one can simply differentiate in the above equation with respect to $x$, then set $x = 0$ to get rid of any other terms not used. Doing this, I get $a_2 = \frac{1}{4}$, $a_3 = \frac{-1}{4}$, $a_4 = \frac{7}{48}$, and $a_5 = \frac{-19}{160}$.

My apologies to Claude!

What we are after is just some $k$'th order Taylor expansion of $y(x)$ about $x=0$. From Taylor's theorem we have that the coefficients $a_n$ in $$f(x) = \sum_{n=0}^\infty a_n x^n$$ satisfy $a_n = \frac{f^{(n)}(0)}{n!}$. We can compute this recursively by taking the $n$'th derivative of the ODE, $y'(x) = -e^{y(x)} - \frac{\sinh(x)}{2}$, to get $$y^{(n+1)}(0) = \left.\frac{d^n}{dx^n}e^{y(x)}\right|_{x=0} - \frac{1}{2}\left.\frac{d^n\sinh(x)}{dx^n}\right|_{x=0}$$ Using the expression for the $n$'th derivative of $e^{y(x)}$ we obtain the recursion

$$a_{n+1} = -\frac{1}{n+1}\left(\sum_{(m_1,\ldots,m_n)} \prod_{j=1}^n\frac{a_j^{m_j}}{m_j!}\right) -\frac{1-(-1)^n}{4(n+1)!}$$

where the sum extends over all $m$-tuples $(m_1,m_2,\ldots,m_n)$ of non-negative integers satisfying $m_1+2m_2+\ldots+nm_n = n$.

The formula above is not very useful when computing it by hand but it can be useful on a computer to automatize the computation of the coefficients. A simple (far from optimal) Mathematica implementation of this is given below and gives the following coefficients for $1\leq i \leq 20$:

$$a_i = \left\{-1,\frac{1}{4},-\frac{1}{4},\frac{7}{48},-\frac{19}{160},\frac{127}{1440},-\frac{185}{2688},\frac{125}{2304},-\frac{84131}{1935360},\frac{255841}{7257600},-\frac{12293681}{425779200},\frac{3263593}{136857600},-\frac{526926271}{26568622080},\frac{2893020049}{174356582400},-\frac{22242690527}{1594117324800},\frac{493179580879}{41845579776000},-\frac{43334418279277}{4335999123456000},\frac{8376548037343}{984980570112000},-\frac{150701255607005953}{20760763803107328000},\frac{465109942159859}{74858523328512000}\right\}$$

nmax = 20;
aa = Table[-1, {i, 1, nmax}];
varlist = Table[Subscript[x, i], {i, 1, nmax}];
Do[
(* Variables *)
var = varlist[[1 ;; i]];

(* Define constraint on tuples *)
constraint = Sum[j varlist[[j]], {j, 1, i}] == i;
Do[constraint = constraint && var[[j]] >= 0, {j, 1, i}];

(* Compute allowed tuples *)
tuples = var /. Solve[constraint, var, Integers];

(* Compute sum over tuples *)
sum = 0;
Do[
curtuple = tuples[[j]];
(* Product over elements in tuples*)
prod = 1;
Do[prod *= 1/curtuple[[k]]! aa[[k]]^curtuple[[k]]];, {k, 1, i}];
sum += prod;
, {j, 1, Length[tuples]}];

(* Store coefficient a(i+1) *)
aa[[i + 1]] = -(sum/(i + 1)) - (1 - (-1)^i)/(4 (i + 1)!);
Print[i + 1, "   ", aa[[i + 1]]];
, {i, 1, nmax - 1}];