- #1
Dustinsfl
- 2,281
- 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:
So to have an array of different dxs, I have written
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:
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
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))
Code:
N = 2 ** np.arange(-4, 10, dtype = np.float64)
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
Here is an image of what I am trying to create:
http://imageshack.us/a/img35/4462/6b67.jpg
Last edited: