Tell me more ×
Computational Science Stack Exchange is a question and answer site for scientists using computers to solve scientific problems. It's 100% free, no registration required.
up vote 4 down vote favorite
1

I am using some finite difference algorithm to solve the problem of a parabolic equation. Reading the Leveque's book on finite differences he suggests to test convergence of the method by considering the ratio of differences between solutions, that is $$\frac{u_h-u_{h/2}}{u_{h/2}-u_{h/4}}$$. I understand that it is valid at any given point on the grid but can I estimate in the same way the convergence of the whole solution which is the vector with a dimension $N$, the number of grid points? My goal is to show convergence of the method where the error $e$ is defined as a difference between numerical solution and the function projected on the grid. That is I want to find $p$ s.t. $||e||=Ch^p$. Let me know if the approach is still valid, thus I would find the ratio not the differences but the ration of those norms, i.e. $$\frac{||u_h-u_{h/2}||}{||u_{h/2}-u_{h/4}||}$$

share|improve this question
he did not say not to use it, he said not to use a method where the reference solution is a very fine grid solution. However, you made a valid point that dimensions don't match, that's right. Could you please suggest me a name of a method/book or anything where I can find how to estimate the convergence of the total error, not at one point error? I don't know the exact solution, I can find only successive numerical solutions. – Kamil Nov 6 '12 at 3:17
basically, I have a theoretical estimate for the method $||e||_h=\sqrt{h\sum_{i=1}^m e_i^2}=O(h^2)$, thus I would expect the error to decrease by 4 every time I refine the mesh twice, I just don't know how to measure this "4". Please advise. – Kamil Nov 6 '12 at 3:28
In principle, you could interpolate both $u_h$ and $u_{h/2}$ on a finer grid (say, $u_{h/4}$) and then compute the difference. However, the recommendation not to do this still stand. – Christian Clason Nov 7 '12 at 13:48
Christian, why interpolate on a finer grid, could not you compare on $h/2$ mesh? Also, when you interpolate,say linearly, you introduce locally second and globally first order error(unless you use a higher degree polynomial)...Also, can you pleae explain why not to use it? I see that very often in thesis etc., where one takes a point on the grid and computes the ratio between differences between solutions. – Kamil Nov 8 '12 at 1:03
The point is precisely to minimize the influence introduced by interpolation. The hope is that by interpolating both solutions, the influence cancels out to some extent (as opposed to only interpolating the coarser one). Of course you'd use higher order polynomials for interpolation than for your solution. – Christian Clason Nov 9 '12 at 11:09
show 1 more comment

4 Answers

up vote 3 down vote accepted
+50

The usual procedure for verifying convergence rates is to take a problem where you know the exact solution analytically. The easiest way to do this is to start from the solution: Pick a function that satisfies the boundary data and then plug it into the differential equation to get the corresponding right-hand side. (This is sometimes called the method of manufactured solution.)

A (very simple) example: Consider the one-dimensional parabolic problem $$ u_t - u_{xx} = f $$ with boundary conditions $u(t,0) = u(t,1) = 0$ and initial conditions $u(0,x) = 0$. Then you could choose, say, $$ u(t,x) = \sin(\pi x)t $$ and compute $$ f(t,x) = u_t(t,x) - u_{xx}(t,x) = (1+\pi^2t)\sin(\pi x)$$ and use this as the right hand side in your finite difference code. This is a bit more complicated if you need to consider, e.g., jumping coefficients, but it can be done by choosing the coefficients together with the right-hand side.

share|improve this answer
This is helpful, thanks! However, I want to test on the rough initial data, because the whole point I am testing how robust the method is. Even though right after time 0 the solution might be infinitely smooth, but I have to start with a solution that has initial data continuous but not differentiable, or even discontinuous....How to I consider this case? Shall I smooth my initial data at the non-differentiable points making sure the radius of smoothness is still smaller than grid size? Just a guess... – Kamil Nov 8 '12 at 1:10
That is of course more difficult, and takes a bit of fiddling. The details depend heavily on the equation you're trying to solve, but you could possibly take your nonsmooth initial data, multiply it with a sufficiently quickly decaying exponential function (in time) and add it to a smooth function with zero initial data. – Christian Clason Nov 9 '12 at 11:25
Christian, in my case it is a heat equation and initial data is a "hokey stick"(the value is $0$ for $x<C$ and a straight line with a slope $k=1$ after), I like this idea with adding the rough piece of data by multiplying say by $e^{-500t}$, I see the impact of that initial data decays rapidly, however, how to manufacture this solution that provides such an initial data? That is as you suggested, I take $u(t,x)=sin(\pi x)t+\max(x-C,0)*e^{-500t}$? Now that I want to calculate $f$ I need to find the derivatives and I can't at $x=C$, how do I deal with it? Please suggest. – Kamil Nov 11 '12 at 4:30
I'm sorry, I should have read your first comment more carefully. If you need to have discontinuous data (as opposed to not infinitely often differentiable), this approach breaks down since you need to consider weak solutions. In this case, it's better to try to construct an exact solution. Maybe this paper helps. – Christian Clason Nov 11 '12 at 12:37
the data is continuous, but not differentiable at one point. However, due to diffusion the solution is infinitely differentiable everywhere inside of the domain, that is it becomes one right after time zero. So inside I am perfectly fine do differentiate it. I have the method, I just want to test it on the convergence now and want to use MMS as I find it promising. At least I can take a smooth function that is in the limit produces the hockey stick... – Kamil Nov 11 '12 at 14:26
show 3 more comments
up vote 3 down vote

Yes, what you're proposing is certainly a valid approach.

I will attempt to also answer what I think you have asked in the comments. If this is the answer you're looking for, then we should clarify the question.

The error is of the form $$\|e_h\| = Ch^p + {\mathcal O}(h^{p+1})$$ where, importantly, $C$ is independent of $h$. Note that if we take the log of the above, we get $$\log(\|e_h\|) \approx \log(C) + p\log(h),$$ so if you plot error versus $h$ on a log-log plot, you should see a straight line with slope $p$.

You can compute some values of $e_h$ for different $h$ using your method. Given two pairs $(h_1, e_{h_1})$ and $(h_2,e_{h_2})$, you can approximate $p$ as follows. Note that $$\frac{\|e_{h_1}\|}{\|e_{h_2}\|} \approx \left(\frac{h_1}{h_2}\right)^p$$ Thus $$p \approx \log\left(\frac{\|e_{h_1}\|}{\|e_{h_2}\|}\right)/\log\left(\frac{h_1}{h_2}\right).$$

share|improve this answer
In the form of the error you are already assuming that the error will be zero on an infinitely fine grid. But this is something that needs proving. A more general form is $\| u_h -u_0\| = C h^p$ where $u_0$ is the value you get for $h\downarrow0$ which is not necessarily $u_\mathrm{exact}$. Now you need at least 3 grid to determine $u_0$, $C$ and $p$. See for example Eca & Hoekstra – chris Nov 9 '12 at 21:10
assume I don't have an opportunity to compute a solution on a very fine grid. Can I at least estimate the error at one point by comparing the ratio between differences in solutions? I see from comments that this is not a good practice but this is what I got from the book and wonder if this is a valid way of estimating the convergence at least at one point. – Kamil Nov 10 '12 at 19:53
Medan : Yes. @chris : You can certainly choose to take a different definition of the word "error", but I think it makes sense to use the one that is implicit in the question. – David Ketcheson Nov 11 '12 at 11:10
@David : if the answer to the question is code verification with the method of manufactured solutions I think it makes sense to use $\|e_h\| = e_0 + Ch^p$ and show that the code gives $e_0$ zero and the expected $p$. Assume the code has a flaw so that $e_0$ is small but not negligible. If you use two relatively coarse grids it could be that $e_0 << \|e_h\|$ and you will not notice the flaw. – chris Nov 11 '12 at 15:54
up vote 1 down vote

The book by Patrick Roach is probably a good reference:

Roache, Patrick J. "Verification and validation in computational science and engineering", Hermosa publishers, 1998.

share|improve this answer
up vote -2 down vote

Isn't it possible to use the root mean square differences between the $u$ values between iterations. The RMS values computed this way should decrease rapidly during the iterations. You could terminate your iterations if the value drops below a threshold.

share|improve this answer
This is not only off-topic, since the question is not about any iterative method, but also a very bad idea in general, since this can mistake stagnation (no progress towards the solution) for convergence. – Christian Clason Nov 7 '12 at 16:17

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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