Solve BVP $y''-3xy'-y=x$, $y(0)=y'(1)=0$ by power series, got differences to numerical solution

Question

I was trying to solve the ODE $y''-3xy'-y=x$ by series. Well, we have regular points all over it so I thought it would be a nice shot, but I am having some problems finding the solutions.

By proposing a solution of the form $y(x)=\sum_{n=0}^{\infty}a_{n}x^{n}$, we will have $y'(x)=\sum_{n=1}^{\infty}a_{n}(n)x^{n-1}$, and $y''(x)=\sum_{n=2}^{\infty}a_{n}(n)(n-1)x^{n-2}$. Substituting in the ODE: \begin{equation} \sum_{n=2}^{\infty}a_{n}(n)(n-1)x^{n-2}- 3x\sum_{n=1}^{\infty}a_{n}(n)x^{n-1}- \sum_{n=0}^{\infty}a_{n}x^{n} =x \end{equation} \begin{equation} \sum_{n=2}^{\infty}a_{n}(n)(n-1)x^{n-2}- \sum_{n=1}^{\infty}3a_{n}(n)x^{n}- \sum_{n=0}^{\infty}a_{n}x^{n} =x \end{equation} \begin{equation} \sum_{n=2}^{\infty}a_{n}(n)(n-1)x^{n-2}- \sum_{n=3}^{\infty}3a_{n-2}(n-2)x^{n-2}- \sum_{n=2}^{\infty}a_{n-2}x^{n-2} =x \end{equation} After, I arrange the equation in a more direct form: \begin{equation} (2a_{2}-a_{0})+ [(3!)a_{3}-4a_{1}]x+ \sum_{n=4}^{\infty}[a_{n}(n)(n-1)x^{n-2}-3a_{n-2}(n-2)x^{n-2}-a_{n-2}]x^{n-2} =x \end{equation} And I find the resulting system: \begin{alignat*}{2} a_{2} & {}={} & \frac{a_{0}}{2} \\ a_{3} & {}={} & \frac{(1+4a_{1})}{3!} \\ a_{n} & {}={} & \frac{3n-5}{n(n-1)}a_{n-2} \end{alignat*} Separating for odd, and even terms:

For even $k\geq2$ \begin{equation} a_{2k}=\frac{(6k-5)}{(2k)(2k-1)}a_{2k-2}=\frac{\Pi_{j=2}^{k}(6j-5)}{(2k)!}a_{0} \end{equation} For odd $k\geq2$: \begin{equation} a_{2k+1}=\frac{(6k-2)}{(2k+1)(2k)}a_{2k-1}=\frac{\Pi_{j=2}^{k}(6j-2)}{(2k+1)!}[1+4a_{1}] \end{equation} And then comes the trick. I want to solve it for the conditions: $y(0)=0$ and $y'(1)=0$. The former says all my even terms are zero (Yey!), so I can abandon them, but the latter says I need to satisfy this relation $\sum_{k=0}^{\infty}a_{2k+1}(2k+1)=0$. Well, in this case I need to have $-0.25<a_{1}<0$, otherwise I couldn't make the sum zero. \begin{equation} \sum_{k=0}^{\infty}a_{2k+1}(2k+1)=a_{1}+3a_{3}+\sum_{k=2}^{\infty}a_{2k+1}(2k+1)=0 \end{equation} \begin{equation} a_{1}+\frac{(1+4a_{1})}{2!}+(1+4a_{1})\sum_{k=2}^{\infty}\frac{\Pi_{j=2}^{k}(6j-2)}{(2k)!}=0 \end{equation} Making $S=\sum_{k=2}^{\infty}\frac{\Pi_{j=2}^{k}(6j-2)}{(2k)!}$, I get that $a_{1}=-\frac{\frac{1}{2}+S}{3+4S}$.

Apparently, I don't know if this is right because I am also solving this same problem using finite differences, and the results are really different one from the other as you can see on the image.

Answer 1

Your formulas all check out. Let's try to replicate the numerical results. Simplify by shifting the inhomogeneity to the first coefficient, set $c=4a_1+1$, then $a_1=\frac{c-1}4$ and the series becomes $$ y(x)=\frac{c-1}4x+c\sum_{k=1}^\infty\frac{\prod_{j=2}^{k}(6j-2)}{(2k+1)!}x^{2k+1} \\ y'(x)=\frac{c-1}4+c\sum_{k=1}^\infty\frac{\prod_{j=2}^{k}(6j-2)}{(2k)!}x^{2k} $$ First evaluate an approximation for $S=\sum_{k=0}^\infty\frac{\prod_{j=2}^{k}(6j-2)}{(2k)!}$ so that for $y'(1)=0$ one needs $$ 0=\frac{c-1}4+cS\iff c=\frac1{1+4S} $$

k=1
term_k = 1/2
S_k = term_k
while S_k+term_k!=S_k :
    k+=1
    term_k *= (3*k-1)/(k*(2*k-1))
    S_k += term_k
c = 1/(1+4*S_k)
n=k
x = np.linspace(0,1.2,150)
term_k=x**3/6
y = term_k
for k in range(2,n):
    term_k *= x**2*(3*k-1)/(k*(2*k+1))
    y+=term_k
y = (c-1)/4*x+c*y
plt.plot(x,y, lw=3, label="series")

For the numerical solution, use the second order approximations of the derivatives and use for expediency the general solver fsolve. Encode the equations and boundary conditions, the right via a ghost cell, and then plot everything

from scipy.optimize import fsolve

for N in [5,10]:
    x = np.arange(N+2)/N
    dx = 1/N
    def F(y):
        y = np.concatenate([[0],y]) # force BC at x=0
        diffeq = (y[2:]-2*y[1:-1]+y[:-2])-3*x[1:-1]*dx/2*(y[2:]-y[:-2]) - dx**2*(y[1:-1] + x[1:-1])
        bc1 = y[-1]-y[-3]
        return [bc1, *diffeq]

    y = fsolve(F,0*x[1:])
    plt.plot(x[1:],y,'-x', lw=0.5,label=f"N={N}")
plt.legend(); plt.grid(); plt.show()

The resulting plot shows that the numerical solution is as close as can be expected to the series solution

As your series solution plot is the same as above, you must have some error in establishing the system for the numerical solution.