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))
 
Suppose ,instead of the usual x,y coordinate system with an I basis vector along the x -axis and a corresponding j basis vector along the y-axis we instead have a different pair of basis vectors ,call them e and f along their respective axes. I have seen that this is an important subject in maths My question is what physical applications does such a model apply to? I am asking here because I have devoted quite a lot of time in the past to understanding convectors and the dual...
Insights auto threads is broken atm, so I'm manually creating these for new Insight articles. In Dirac’s Principles of Quantum Mechanics published in 1930 he introduced a “convenient notation” he referred to as a “delta function” which he treated as a continuum analog to the discrete Kronecker delta. The Kronecker delta is simply the indexed components of the identity operator in matrix algebra Source: https://www.physicsforums.com/insights/what-exactly-is-diracs-delta-function/ by...

Similar threads

Back
Top