I want to solve a Sturm-Liouville problem in 1D, i.e., \begin{align} [p(x)\ u'(x)]'+q(x)u(x) = f(x) \end{align} with boundary conditions \begin{align} u(0)=a \hspace{1cm} u'(0)=b \end{align} How do I implement these conditions?

What I did:

First I multiply both sides by a test function $v\in\lbrace c\in H^1\big\vert v(0)=0, v'(0)=0 \rbrace$ and integrate over the computational domain, i.e. from 0 to 1 in my case. \begin{align} -\int\limits_0^1\left(p(x)u'(x)\right)'v(x)dx +\int\limits_0^1q(x)u(x)v(x)dx = \int\limits_0^1f (x)v(x)dx \end{align} The first term can be transformed by partial integration. \begin{align} -p(1)u'(1)v(1) + \underset{=0}{\underbrace{p(0)u'(0)v(0)}} +\int\limits_0^1p(x)u'(x)v'(x)dx \end{align} I want to use hermite cubic splines denoted by $\phi_k^l$ with the property \begin{align} \left(\phi_k^l\right)^{(m)}(x_j)=\delta_{kj}\delta_{lm} \end{align} where $(m)$ denotes the 0th or 1st derivative. Inserting this in the weak formulation above we finally arrive at the matrix equation \begin{align} \nonumber -p(1)u^1_{N +1}&\delta_{0l}\delta_{N+1j} +a\int\limits_0^1 dx p(x)\left[\phi_0^0\,'(x)\phi_j^l\,'(x) + q(x)\phi_0^0(x)\phi_j^l(x)\right]\\ & +b\int\limits_0^1 dx\, p(x)\left[\phi_0^1\,'(x)\phi_j^l\,'(x) + q(x)\phi_0^1(x)\phi_j^l(x)\right]\\ &+\sum\limits_{i=1}^{N+1}\sum\limits_{k=0}^1u_i^k \int\limits_0^1 dx\, p(x)\left[\phi_i^k\,'(x)\phi_j^l\,'(x) + q(x)\phi_i^k(x)\phi_j^l(x)\right]\\ =&\int\limits_0^1f (x)\phi_j^l (x)dx \hspace{2cm}j=1,\dots,N+1\ \ l=0,1 \end{align} Is this derivation correct up to here?

The matrix equation is set up by first mapping every element to the standard intervall $[-1,1]$ and then using Gauss Integration.

Afterwards the global matrix is assembled. Finally boundary conditions have to be implemented according to the above equation. This means we delete the lines and columns corresponding to the coeffiicients of $\phi_0^0$ and $\phi_0^1$. Subtracting the second and third term (the ones with a and b in front) from the left hand side of the above equation on the right hand side and subtraction $p(1)$ at the corresponding position of the matrix we arrive at the final matrix equation. This is visualized in the figure Implementation of boundary conditions.

But using the extremely simple case of $p(x)=q(x)=1$, $f(x)=0$ and using only one domain we can find an analytic solution. The comparison between the analytic and the numerical solution is shown in figure comparision of solutions. Thus there must be a mistake in the above implementation process (other boundary conditions provided at both ends of the domain give correct solutions.)

For an increasing number of domains the solution does also not converge. Therefore something is obviously wrong but I do not have any idea.

Can you please help me with this?

  • $\begingroup$ I believe you have a small error in your figure showing implementation of the boundary conditions-- the column multiplied by b should be the one to the left. Otherwise, I agree with what you have done. I did an implementation of this and I am seeing a converged solution with a large number of elements (say 100) but to the wrong solution! I don't yet understand why but would be interested to know if you see the same? $\endgroup$ Nov 1 '15 at 17:24
  • $\begingroup$ You are right, there is a typo in the figure showing the implementation of the boundary conditions, it should be $K_{32}$ and $K_{42}$. For a larger number of elements I am also seeing convergence to the wrong solution. The numerical solution behaves as $u(1) \rightarrow 0$. As this is my first larger numerical project I have no expirience thus I am not sure if that helps, but I also used Spectral Element Methods (Chebyshev Polynomials in several subintervals) for this problem and again see the tendency towards $u(1) \rightarrow 0$. $\endgroup$
    – NoIdea
    Nov 3 '15 at 11:22
  • $\begingroup$ I am not sure what you mean by "The numerical solution behaves as u(1)→0"? There is no BC at u(1) so how can you make the solution approach zero there? $\endgroup$ Nov 3 '15 at 15:33
  • $\begingroup$ For a large number of elements the obtained solution (not the correct analytic one) seems to be 0 at the endpoint $x=1$, i.e. the coefficient for the last basis function seems to be tend to 0. That's what I meant by $u(1)\rightarrow 0$. So I do not make the solution to be zero there but the numerical outcome turns to be zero. $\endgroup$
    – NoIdea
    Nov 3 '15 at 15:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.