Plotting the error vs the change in 1/dx

[Dustinsfl]

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:

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.0for 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

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:

  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:
http://imageshack.us/a/img35/4462/6b67.jpg

 

[Dustinsfl]

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

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))