# Plotting the error vs the change in 1/dx

#### dwsmith

##### Well-known member
In order to plot the error vs the change in 1/dx, I need to integrate the the solution with different a dx and take the max error. However, I don't know how I can set this up in python or matlab.

In python, my code for just one dx is:
Code:
import numpy as np
L = 80.0
N = 512
dt = 0.0002
tmax = 10
nmax = int(np.floor(tmax / dt))
dx = L / N
x = np.arange(-L / 2.0, L / 2.0 - dx, dx)
k = np.hstack((np.arange(0, N / 2.0 - 1.0),
np.arange(-N / 2.0, 0))).T * 2.0 * np.pi / L
k1 = 1j * k
k3 = (1j * k) ** 3
u = 2 * (2 / (np.exp(x + 20.0) + np.exp(-x - 20.0))) ** 2
udata = u
tdata = 0.0

for nn in range(1, nmax + 1):
du1 = (-np.fft.ifft(k3 * np.fft.fft(u)) -
3 * np.fft.ifft(k1 * np.fft.fft(u ** 2)))
v = u + 0.5 * du1 * dt
du2 = (-np.fft.ifft(k3 * np.fft.fft(v)) -
3 * np.fft.ifft(k1 * np.fft.fft(v ** 2)))
v = u + 0.5 * du2 * dt
du3 = (-np.fft.ifft(k3 * np.fft.fft(v)) -
3 * np.fft.ifft(k1 * np.fft.fft(v ** 2)))
v = u + du3 * dt
du4 = (-np.fft.ifft(k3 * np.fft.fft(v)) -
3 * np.fft.ifft(k1 * np.fft.fft(v ** 2)))
u = u + (du1 + 2.0 * du2 + 2.0 * du3 + du4) * dt / 6.0
if np.mod(nn, np.ceil(nmax / 20.0)) == 0:
udata = np.vstack((udata, u))
tdata = np.vstack((tdata, nn * dt))
So to have an array of different dxs, I have written
Code:
N = 2 ** np.arange(-4, 10, dtype = np.float64)
but then how do I adjust the rest of code? Or is there a better way to be able to obtain the error plot which is of the form $$\exp(-c\cdot dx)$$?

If python doesn't suit you, I have matlab code too:
Code:
  L = 80;
N = 512;
dt = 0.0002;
tmax = 10;
nmax = round(tmax/dt);
dx = L/N;
x = (-L/2:dx:L/2-dx)';
k = [0:N/2-1 -N/2:-1]'*2*pi/L;
k1 = 1i*k;
k2 = (1i*k).^2;
k3 = (1i*k).^3;
u = 2*sech(x + 20).^2;
udata = u;
tdata = 0;

% integration begins
for nn = 1:nmax
du1 = -ifft(k3.*fft(u)) - 3*ifft(k1.*fft(u.^2));
v = u + 0.5*du1*dt;
du2 = -ifft(k3.*fft(v)) - 3*ifft(k1.*fft(v.^2));
v = u + 0.5*du2*dt;
du3 = -ifft(k3.*fft(v)) - 3*ifft(k1.*fft(v.^2));
v = u + du3*dt;
du4 = -ifft(k3.*fft(v)) - 3*ifft(k1.*fft(v.^2));
u = u + (du1 + 2*du2 + 2*du3 + du4)*dt/6;
if mod(nn, round(nmax/20)) == 0
udata = [udata u];
tdata = [tdata nn*dt];
end
end
% integration ends
I have even less of a clue on what do for the matlab code.

Here is an image of what I am trying to create:

Last edited:

#### dwsmith

##### Well-known member
I have progressed some with the code in Python but its having a shape problem. Hopefully some can help square it way.

Code:
import numpy as np
import pylab

L = 80.0
dt = 0.0002
tmax = 10
nmax = int(np.floor(tmax / dt))
deltax = []
error = []

for N in [80., 100., 120., 140., 160., 180., 200., 220., 240., 260., 280., 300.,
320., 340., 360.]:
dx = L / N
deltax.append(dx)
x = np.arange(-L / 2.0, L / 2.0 - dx, dx)
k = np.hstack((np.arange(0, N / 2.0 - 1.0),
np.arange(-N / 2.0, 0))).T * 2.0 * np.pi / L
k1 = 1j * k
k3 = (1j * k) ** 3
u = (2 * (2 / (np.exp(x + 20.0) + np.exp(-x - 20.0))) ** 2)
udata = u
tdata = 0.0
#  integration
for nn in range(1, nmax + 1):
du1 = (-np.fft.ifft(k3 * np.fft.fft(u)) -
3 * np.fft.ifft(k1 * np.fft.fft(u ** 2)))
v = u + 0.5 * du1 * dt
du2 = (-np.fft.ifft(k3 * np.fft.fft(v)) -
3 * np.fft.ifft(k1 * np.fft.fft(v ** 2)))
v = u + 0.5 * du2 * dt
du3 = (-np.fft.ifft(k3 * np.fft.fft(v)) -
3 * np.fft.ifft(k1 * np.fft.fft(v ** 2)))
v = u + du3 * dt
du4 = (-np.fft.ifft(k3 * np.fft.fft(v)) -
3 * np.fft.ifft(k1 * np.fft.fft(v ** 2)))
u = u + (du1 + 2.0 * du2 + 2.0 * du3 + du4) * dt / 6.0
#    error.append(max(abs(udata[:,-1] - 2. *
#                     (2. / (np.exp(x - 20) + np.exp(-x - 60))))))
if np.mod(nn, np.ceil(nmax / 20.0)) == 0:
udata = np.vstack((udata, u))
tdata = np.vstack((tdata, nn * dt))