Discrete solutions of boundary value problem with parameter.

Solution 1:

This is a linear Sturm–Liouville eigenvalue problem (if $f'(0) = 1$ condition is dropped, it is only needed to fix the unique solution). $$ -\frac{d}{dx}\left[p(x) \frac{df}{dx}\right] + q(x) f(x) = \lambda w(x) f(x) $$

First we need to scale the equation by $C(x)$ to form the $-\frac{d}{dx}\left[p(x) \frac{df}{dx}\right]$ term. $$ C(x) (x-3)^7 ((x^2-4 x+3) f''(x)+2 x f'(x)) = -\frac{d}{dx}\left[p(x) \frac{df}{dx}\right] $$ Eliminating $C(x)$ gives an equation for $p(x)$ $$ \frac{2x}{x^2-4 x+3} = \frac{p'(x)}{p(x)} = (\log p(x))' $$ By integrating we obtain $$ p(x) = \frac{(x-3)^3}{x-1} $$ and $$ C(x) = \frac{-p(x)}{(x-3)^7 (x^2 - 4x + 3)} = -\frac{1}{(x-1)^2(x-3)^5}. $$

$$-\frac{d}{dx}\left[\frac{(x-3)^3}{x-1} \frac{df}{dx}\right] -\frac{3456}{(x-1)^2(x-3)^5} k (3-2 x)^2 (x-1) x^2 f(x)=0, \quad x \in [0,1]$$ Moving term with $k$ to the right we obtain the Sturm-Liouville form $$ -\frac{d}{dx}\left[\frac{(x-3)^3}{x-1} \frac{df}{dx}\right] =k \cdot \frac{3456}{(x-1)(x-3)^5} (3-2 x)^2 x^2 f(x)=0, \quad x \in [0,1] $$ with parameters $$ p(x) = \frac{(x-3)^3}{x-1}, \quad q(x) = 0, \quad w(x) = \frac{3456 x^2 (3-2 x)^2}{(x-1)(x-3)^5}. $$

The idea of the numerical method is to introduce a regular grid on $x \in [0, 1]$ and approximate the equation with second order as $$ -\frac{p(x_{i+1/2}) (f_{i+1} - f_i) - p(x_{i-1/2}) (f_i - f_{i-1})}{h^2} + q(x_i) f_i = k w(x_i) f_i, \quad i = 1, \dots, N-1\\ f_0 = f_N = 0 $$ where $h = \frac{1}{N}$ and $x_j = jh$.

This results in an (generalized) eigenvalue problem $$ A \mathbf f = k W \mathbf f $$ with tridiagonal matrix $A$ and diagonal matrix $W$, both symmetric and positive definite. The vector $\mathbf f = (f_1, f_2, \dots, f_{N-1})$. Nonzero elements of $A$ and $W$ are $$ A_{i,i+1} = -p(x_{i+1/2})\\ A_{i,i-1} = -p(x_{i-1/2})\\ A_{i,i} = p(x_{i-1/2}) + p(x_{i+1/2}) + h^2 q(x_i)\\ W_{i,i} = h^2 w(x_i). $$ The original eigenvalues $k$ can be approximated as eigenvalues of $W^{-1/2} A W^{-1/2}$.

The following Mathematica code implements the method

p[x_] = (3-x)^3/(1-x);
w[x_] = (3456x^2 (2x-3)^2)/(3-x)^5/(1-x);
q[x_] = 0;

M = 5000;
h = 1/M;
x[j_] = j h;

B = SparseArray[{i_, j_} /; Abs[i - j] <= 1 :> 
    If[i == j, (p[x[i - 1/2]] + p[x[i + 1/2]] + h^2 q[x[i]])/(h^2 w[x[i]]), 
               -p[x[(i + j)/2]]/(h^2 Sqrt[w[x[i]] w[x[j]]])
    ], {M - 1, M - 1}];

Sort@Eigenvalues[N@B,-10]
// {4.68175, 16.7411, 36.1904, 63.0365, 97.2817, 138.927, 187.973,
244.42, 308.269, 379.519}
```