How to build linear system of approximation of solution to boundary value problem

Question

We want to numerically approximate the following problem using finite differences.

$$y''=f(t,y,y') \ \ \ \ \ a \leq t \leq b$$ $$y(a) = \alpha \ \ \ \ \ y(b) = \beta$$

We can divide the interval $[a,b]$ into the subintervals defined by the points

$$t_i = a+ih; \ \ \ h = \frac{b-a}{n}$$

where $t_0 = a$ and $t_n =b$.

The finite differences approximations give

$$y'(t_i) \approx \frac{y_{i+1} - y_{i-1}}{2h}$$ $$y''(t_i) \approx \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2}$$

Therefore, the problem becomes the "approximated" problem

$$y'' - f(t,y,y') = 0 \to \frac{y_{i+1} - 2y_i + y_{i-1}}{h^2} - f(t_i, y_i, y'_i) = 0$$

$$y_{i+1} - 2y_i + y_{i-1} - h^2f \bigg(t_i, y_i, \frac{y_{i+1} - y_{i-1}}{2h} \bigg) = 0$$

I know this is supposed to be a linear system representable with a tridiagonal matrix, but I fail to find out how to construct this system explicitly.

Moreover, what am I supposed to do if the problem has Neumann conditions on one or both of the boundary points? Or if I have a vector-valued problem?

Answer 1

If you have some approximation $y_i$ that is at least qualitatively correct, then a correction $y_i-u_i$ to it can be computed from the linearization of the problem. \begin{multline} 0=\left(\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}-f\left(t_i,y_i,\frac{y_{i+1}-y_{i-1}}{h}\right)\right) \\ -\left(\frac{u_{i+1}-2u_i+u_{i-1}}{h^2}-\frac{\partial}{\partial y}f\left(t_i,y_i,\frac{y_{i+1}-y_{i-1}}{h}\right)u_i-\frac{\partial}{\partial y'}f\left(t_i,y_i,\frac{y_{i+1}-y_{i-1}}{h}\right)\frac{u_{i+1}-u_{i-1}}{h}\right) \end{multline} that is, you get equations $$ (1-hb_i)u_{i+1}-(2+h^2a_i)u_i+(1+hb_i)u_{i-1}=h^2c_i $$ where \begin{align} a_i&=\frac{\partial}{\partial y}f\left(t_i,y_i,\frac{y_{i+1}-y_{i-1}}{h}\right)\\ b_i&=\frac{\partial}{\partial y'}f\left(t_i,y_i,\frac{y_{i+1}-y_{i-1}}{h}\right)\\ c_i&=\frac{y_{i+1}-2y_i+y_{i-1}}{h^2}-f\left(t_i,y_i,\frac{y_{i+1}-y_{i-1}}{h}\right) \end{align}

The solution $u_i$ of that system is then the Newton update in $y_i$, that is, $y_i-u_i$ is the improved value, attached with all the convergence properties and non-convergence possibilities of the Newton method.

---

For Neumann boundary conditions you can either add "exterior" points $(t_{-1},y_{-1})$ resp. $(t_{n+1},y_{n+1})$ and use the central difference quotient in addition to the ODE in those points; or use the forward and backward differentiation formulas of second order.

---

As test case take a non-linear oscillator $$y''+0.1y'+y+y^3=0.5\cos(t), \text{ that is } ~~ f(t,y,y')=\cos(t)-0.1y'-y^3$$ with arbitrary bc $y(0)=1$, $y'(\pi)=-1$. Take $y(t)=\cos(t/2)-1$ as initial approximation.

plots, graphs = 2, 4
f, ax = plt.subplots(plots, 1, figsize=(10,5*plots))

k, R = 0.1, 0.5
t0, tf = 0, pi
y0, Dyf = 1, -1 
n = 20
h = (tf-t0)/n
t = np.arange(n+2)*h
y = -1+2*np.cos(t/2)
Dy  = np.zeros([n+2])
D2y = np.zeros([n+2])
M = np.zeros([n+2,n+2])
C = np.zeros([n+2])
row = -1
for k in range(plots*graphs):
    Dy[1:-1]  = (y[2:]-y[:-2])/(2*h)
    D2y[1:-1] = (y[2:]-2*y[1:-1]+y[:-2])/h**2
    M[0,0] = 1; C[0]=y[0]-y0
    for i in range(1,n+1):
        a, b, c = -1-3*y[i]**2, -k, D2y[i] + k*Dy[i] + y[i]+y[i]**3-R*cos(t[i])
        M[i,i-1:i+2] = [ 1-0.5*h*b, -(2+h**2*a), 1+0.5*h*b ]
        C[i] = h**2*c
    M[n+1,n-1:] = [-1, 0, 1]; C[n+1] = 2*h*(Dy[n]-Dyf)
    u = np.linalg.solve(M,C)
    y = y-u;
    if k % graphs == 0: row = row+1; dia = ax[row]
    dia.plot(t,y,label="iter=%d"%k)
for dia in ax: dia.legend(); dia.grid(); #dia.set_ylim([-5,10])
plt.show()

You can play with the constants in the equation to get cases where the iteration wildly diverges.