0
$\begingroup$

Background

I am working with chapter 2 in LeVeque's book: https://faculty.washington.edu/rjl/fdmbook/

I have build my own solver in Python to solve the 2 point BVP: $$ \epsilon u''+u(u'-1) =0 , \\ u(0)=\alpha, u(1)= \beta $$

I have followed the axact steps as described by the answer by VoB in this post: Non-linear Boundary Value Problem. How to compute the Jacobian? (uniform grid; Newton method solve for $G(U)=0$ wrt U).

For the parameters $[\epsilon, \alpha, \beta ] = [0.01, -1,1.5]$ I get this result

enter image description here

In the 10th run it has converged to what I believe is the true solution. This looks pretty much correct to me. (run1 is the intial guess in Newton steps). My stepsize is $\Delta x = h = 0.005$

MY QUIESTION

In the post it is mentioned that the order of convergence is 2. But how do I compute the errors to determine the order of convergence to be 2?

Improve this question
$\endgroup$
3
  • $\begingroup$ See math.stackexchange.com/questions/1286926/… for a mathematical treatment of the solution approximation. The jump in the first order approximation $$u(x)=-1+x+a\tanh(a(x-x_0)/(2ϵ))$$ is expected to be where $$(-1+x_0)+a = 0 = (0.5+x_0)-a,$$ so that $a=0.75$ and $x_0=1-a=a-0.5=0.25$. $\endgroup$
    – Lutz Lehmann
    Commented Mar 21, 2021 at 20:28
  • $\begingroup$ @LutzLehmann Maybe it is just me who is a bit stupid here, but how I can use this to compute errors and order of convergence? Thanks $\endgroup$
    – k.dkhk
    Commented Mar 22, 2021 at 10:20
  • $\begingroup$ If that were the case, I would have posted it as an answer. It just shows that your solution is qualitatively correct. And could provide an almost correct initialization for the solver. /// I do not know why the existing answer was deleted, the step-halving or doubling strategy is valid. The only additional point is that one needs to use some function space norm, $L^1$ or $L^2$, to even out the larger errors at the rapid transition in the boundary layer. These would dominate a supremum norm, which makes it less suitable for the error analysis. $\endgroup$
    – Lutz Lehmann
    Commented Mar 22, 2021 at 10:32

1 Answer 1

Reset to default
1
$\begingroup$

The way I computed the solution in the linked answer is the classical one: I just took a reference solution (assuming the code was correct, i.e. the numerical solution I found is the right one) with a small enough step size $h$, say $h=10^{-9}$. Then I computed the solution with smaller $h$s and for each one of those $h$ I computed the $|| e_h ||$ in a double logarithmic plot. I think in that case I used the supremum norm, because the choice of the parameters didn't produce boundary layers. As pointed out by Lutz in comments, this norm might not be the best one for that situations.

Be careful, this is telling you that the method is implemented correctly, not that your solution is the right one. If there are typos in your setup, you might be converging to another solution, not to the desired one.

I think you can use the Octave/MatLab code I provided in that answer to test the order of convergence. If you want to stick with Python, just define something like

m_range = [10**(i) for i in range(0,4)] #number of points in each grid

and loop over this. In the body, just compute the solution and check the error norm, appending this to a list, say error_list.

In the end, just do plt.loglog(m_range,error_list) or stuff like that. It's not important that your line fits the dots of the error, as in my answer (there I used a trick to shift it towards the error dots): the important point is that the error dots are parallel to the line with slope $-2$.

Improve this answer
$\endgroup$
4
  • $\begingroup$ Thank you very much! But how exactly do I find a an error. Lets say I compare $h^{-1}$ with $h^{-9}$. Do I find a certain point in $x^*$ that they both share take the difference or what else? $\endgroup$
    – k.dkhk
    Commented Mar 22, 2021 at 15:21
  • 1
    $\begingroup$ Let's say that your reference solution, u_ref has been computed with m_ref points. Now assume that you're inside the loop over m_range, and you computed a solution u with m points. To compute the infty norm of the error, you can do norm(u-u_ref(1:(m_ref-1)/(m-1):m_ref),inf) What happens is that you're comparing the solutions at the correct points. @k.dkhk $\endgroup$
    – VoB
    Commented Mar 22, 2021 at 16:04
  • $\begingroup$ @k.dkhk Did you manage to solve the task? $\endgroup$
    – VoB
    Commented Mar 22, 2021 at 18:45
  • 1
    $\begingroup$ Not Yet. But I will be looking at it again tomorrow this time. I will be happy to let you know how it goes. $\endgroup$
    – k.dkhk
    Commented Mar 22, 2021 at 18:51

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.