Start this weeks assignments by updating your repository! Either press the sync fork button on gihub.com, or pull in the latest changes through
git pull
This should pull in a new file VibFD.py containing a new test and a vibration equation solver. The solver is implemented as a class in order to experiment easily with different numerical schemes, to be implemented in the __call__ method. The class also contains methods for computing the l2-error and the order of convergence.
The original solver described in the first chapter of Finite Difference Computing with PDEs is implemented in the class VibHPL. In addition to this solver I have created three more solver classes VibFD2, VibFD3 and VibFD4. The first two are second order and the last should be fourth order. The assignment is to implement the `__call__` method in these three classes such that the test in test_order will pass. This tests the order of the implemented solvers.
i) For VibFD2 use a central difference scheme
\(\frac{u^{n+1}-2u^n+u^{n-1}}{\Delta t^2} + \omega^2 u^{n}=0, \quad n\in (1, 2, \ldots, N-1)\)
and then use Dirichlet boundary conditions for both sides of the domain: \(u(0)=I\) and \(u(T) = I\), where \(T=N\Delta t\). Dirichlet boundary conditions are set by modifying the first and last rows of the assembled matrix.
ii) For VibFD3 use the same central difference scheme, but modify the boundary conditions to: \(u(0)=I\) and \(u'(T)=0\). Here you only need to modify the last row of the assembled matrix from i).
iii) For VibFD4 use a fourth order central difference scheme
\(\frac{-u^{n+2} + 16u^{n+1}-30u^{n}+16 u^{n-1}-u^{n-2}}{12\Delta t^2} + \omega^2 u^{n}=0, \quad n\in (2, \ldots, N-2) \)
For \(n=1\) and \(n=N-1\) use a skewed scheme:
\(\frac{u^{n+4}-6 u^{n+3}+14 u^{n+2}- 4u^{n+1}-15 u^{n} + 10 u^{n-1}}{12 \Delta t^2}\)
\(\frac{10 u^{n+1} - 15 u^n - 4 u^{n-1} + 14 u^{n-2} - 6 u^{n-3}+u^{n-4}}{12 \Delta t^2}\)
and set Dirichlet boundary conditions exactly like in i) for rows 0 and N.
Note that it is a little bit more difficult to prove fourth order convergence than second. This is because the error disappears fast, and also because even though the error in the finite difference method is leading fourth order, this does not mean that the fifth order is negligible. And the fifth order error is closer to the fourth than the third order is to the second. So you may need to use higher tolerance (tol) for the test_order function.
iv) Finally, we would like to extend the solvers to work for any analytical (manufactured) solution. To this end we need to solve the vibration equation with a non-zero right hand side
\(u''+\omega^2 u = f\)
where \(f(t)\) will depend on the chosen solution \(u(t)\). Modify the solver VibFD2 such that the modified vibration equation is solved and verify the implementation. It should still be second order accurate. Test using for the method ue both \(u(t)=t^4\) and \(u(t)=\exp(\sin(t))\). Note: do not use a too big domain and set the boundary conditions using the exact \(u(0)\) and \(u(T)\) . This way T does not need to be n*pi/w, which is required for the exact cosine solution.
For this equation you should use sympy, both to set the exact solution and to compute \(f(t)\). The numerical scheme is
\(\frac{u^{n+1}-2u^n+u^{n-1}}{\Delta t^2} + \omega^2 u^{n}=f^n, \quad n\in (1, 2, \ldots, N-1)\)
Note: This assignment is exactly like i), except that you need to include \(f^n\) in the assembled right hand side vector.