1

I wrote a program in Python to solve the 1D wave equation by treating it as a first-order in time linear system, where the spatial variables are discretized (using 2nd order methods) and I integrate in time using a 4th-order Runge-Kutta.

My problem is: when I compute the order of convergence at each time-step, it is roughly equal to 2 for small times, and then gets worse for larger times. Here is a plot (on the y-axis is the L^2 norm of the order of convergence): Order of convergence

In the above I used a spatial grid in the interval [0, 1] with grid interval 0.5, and 10000 time iterations for Runge-Kutta. I computed the solution at the final time T = 0.8 s and thus the time interval was T/10000. I considered periodic boundary conditions for the solutions.

My code to compute the order of convergence is the following:

def orderOfConvergence(u, u2, u3, length, f):
    num = 0
    den = 0
    for i in range(0, length):
        den += (u[i] - u2[int(i/f)])**2
        num += (u2[int(i/f)] - u3[int(i/(f**2))])**2
    if den > 0 and num > 0:
        return 0.5*log(num/den, f)
    else:
        return 0

where u, u2 and u3 are the three solutions computed with finer and finer grids. In my code, f=0.5, since halved the grid interval each time.

What might be the cause for such a behaviour?

moonknight
  • 334
  • 1
  • 7
  • 1
    Your time step might be too small. For RK4 in problems with mediums sized variables and constants, the optimal step size is around 1e-3. For smaller step sizes the accumulation of floating point errors dominates. – Lutz Lehmann Jan 16 '22 at 15:08
  • @LutzLehmann thank you! Indeed I tried with bigger time steps and the order of convergence remains close to 2. The thing is: I could only try this when the final time is smaller than 0.3 s . If I increase the value of the final time I actually need to shrink the time step or the solution will not converge anymore. So maybe the problem was what you suggested, together with the fact that I was considering too large final times? – moonknight Jan 16 '22 at 15:39
  • 1
    All PDE, if considered as ODE with space-differential operators on the right side, are stiff. The discretization inherits that to some degree, unfortunately increasingly with smaller space-step-sizes. So it is to be expected that with explicit methods, the high-frequency components diverge rapidly. To avoid that, one has to use implicit methods or exponential methods. – Lutz Lehmann Jan 16 '22 at 16:09
  • The latter is operator-splitting methods using the matrix exponential on the linear contents. In many cases this can be implemented relatively efficiently using FFT on the space discretization. See https://stackoverflow.com/questions/29803342 and links there for an example. – Lutz Lehmann Jan 16 '22 at 16:09

0 Answers0