Accuracy

[trixxx3]

Im trying to solve a nonlinear ODE with this package and encountered a problem with th accuracy parameter. Following simple example should demonstrate:

import findiff as fd
import numpy as np
import matplotlib.pyplot as plt

plt.figure('low res')
x = np.linspace(-np.pi*5,np.pi*5,100)
y = np.sin(x)
plt.plot(x,y,'r-')
plt.plot(x,y,'ro')
dx = x[1]-x[0]

#init differentiateoperators
d_dx = fd.FinDiff(0,dx,1,acc=2)
dw_dx = d_dx(y)
plt.plot(x, dw_dx, label='low res, low acc')
d_dx = fd.FinDiff(0,dx,1,acc=60)
dw_dx = d_dx(y)
plt.plot(x, dw_dx, label='low res, high acc')
plt.ylim((-1.1,2))
plt.legend()

plt.figure('high res')
x = np.linspace(-np.pi*5,np.pi*5,2500)
y = np.sin(x)
plt.plot(x,y,'r-')
plt.plot(x,y,'ro')
dx = x[1]-x[0]

#init differentiateoperators
d_dx = fd.FinDiff(0,dx,1,acc=2)
dw_dx = d_dx(y)
plt.plot(x, dw_dx, label='high res, low acc')
d_dx = fd.FinDiff(0,dx,1,acc=60)
dw_dx = d_dx(y)
plt.plot(x, dw_dx, label='high res, high acc')
plt.ylim((-1.1,2))
plt.legend()

Of Course a high accuracy at low resolution should lead to an unaccurate results since the stencil is to long for a good result. but i am stunned how inappropiate the result for the high accuracy high resolution is in the central part.

Even though the forward and backward schemes are pretty good, the central scheme is still far off obviously visible. Since the central scheme is mor sccurate as lomng as i understood FD, why isnt the result as prescise as it is at the border?

Do I missunderstand something? Or do i use it wrongly?

I hope i can get help. But it is a really nice and convienient package for Simulations ;-) All the Best trixxx3

[maroba]

Hi,

the problem is that you are trying to use an extremely high accuracy order (60) on an equidistant grid. The finite difference method is inherently unstable in that case. The reason is the so-called Runge phenomenon, which also appears for Lagrange interpolators of high order. To avoid that, you would have to use a non-equidistant grid. With finite difference methods, you would normally want to find the sweet spot between performance, accuracy and stability. Personally, in practice, I wouldn't go higher than order 6. If you need higher accuracy orders, for instance because you can only afford to calculate on a coarse grid, then you should consider spectral methods. If you have periodic boundary conditions, then Fourier methods would be the way to go. If you have a non-periodic case, Chebyshev methods are a great spectral method. The cost for that is to use a non-equidistant grid, though. For a great practical introduction to spectral methods, I can highly recommend "Spectral Methods with MATLAB" by Lloyd Trefethen.

Best regards, Matthias

[trixxx3]

Thanks for the fast reply Matthias,

okay than it is a general issue i cannot avoid. The gererall thing is, that i try to reproduce a papers outcome on an equation like dw/dt = A * d^4w/dx^4 + B * d^2w/dx^2 + C * w using FowardTimeCentralSpace method and i cannot set dx and dy to their default values since it runs unstable in my case thats is why i tried higher accuracy to reduce the errors.

Lets see if i can handle Spectral Methods. I tried to avoid it since my time in understanding a new technique is limited.

Best regards Lukas