MHB Plotting the error vs the change in 1/dx

AI Thread Summary
To plot the error versus the change in 1/dx, the integration solution needs to be adjusted for multiple values of dx in both Python and MATLAB. In Python, the user has defined an array of N values to calculate different dx values, but requires guidance on how to modify the integration loop to accommodate these varying dx values and capture the maximum error. The MATLAB code faces similar challenges, needing adjustments to handle multiple dx values while maintaining the integration process. The goal is to generate an error plot that follows the expected exponential decay pattern with respect to dx. Further assistance is sought to resolve shape issues in the Python code and to clarify the MATLAB implementation.
Dustinsfl
Messages
2,217
Reaction score
5
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.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
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:
http://imageshack.us/a/img35/4462/6b67.jpg
 
Last edited:
Mathematics news on Phys.org
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))
 
Seemingly by some mathematical coincidence, a hexagon of sides 2,2,7,7, 11, and 11 can be inscribed in a circle of radius 7. The other day I saw a math problem on line, which they said came from a Polish Olympiad, where you compute the length x of the 3rd side which is the same as the radius, so that the sides of length 2,x, and 11 are inscribed on the arc of a semi-circle. The law of cosines applied twice gives the answer for x of exactly 7, but the arithmetic is so complex that the...

Similar threads

Back
Top