Due Tuesday March 1, 2016. Do the homework by modifying this Julia notebook file. Your complete homework should have Julia code in code cells followed by plots that illustrate correct working behavior of the code. Submit the homework as a Julia notebook file to me via email or Canvas.
Write a numerical simulation algorithm for the Kuramoto-Sivashinksy equation $u_t = -u_{xx} - u_{xxxx} - u u_x$ on a periodic spatial domain $x \in [0,L]$ with periodic boundary conditions, using Fourier discretization in space and Crank-Nicolson/Adams-Bashforth semi-inplicit temporal discretization. Use the initial condition $u(x,0) = \cos(ax) + 0.01 \cos(2ax)$ where $a=2\pi/L$, parameter values $r = 0.2, L_x = 100, N_x = 256$, and $\Delta t = 1/16$, and integrate from $t=0$ to $t=100$.
Make a colorplot of $u(x,t)$ in the $x,t$ plane, using every x gridpoint but plotting t at a larger interval than $\Delta t$, perhaps at intervals of 1/2 or 1. It should look just like the plot below.
In [1]:
#Parameter definitions
Lx=100
t=100
Nx=256
dt=1/16
dx=float(Lx/Nx)
x=collect(0:dx:Lx-dx)
T=collect(0:dt:t-dt)
Nt=round(Int,t./dt)
Out[1]:
In [3]:
a=2.5*pi/Lx
u=cos(a*x)+0.01*cos(2*a*x);
In [4]:
function cnabks(u,Lx,T,Nx,dt)
kx=vcat(0:Nx/2-1,0,-Nx/2+1:-1) #integer wavenumbers
alpha=2*pi*kx/Lx #non-dimensionalized wavenumbers
#operators in Fourier basis
D=1im*alpha #D_x(e^{i alpha x}=i alpha e^{i alpha x}
L=-D.^2-D.^4 #Lu=-u_{xx}-u_{xxxx}
G=-.5*D #(1/2)D_x u^2= u u_x
hatu=zeros(Complex{Float64},Nx,Nt)
hatNu=zeros(Complex{Float64},Nx,Nt)
usol=zeros(Complex{Float64},Nx,Nt)
hatu0=fft(u)
hatu2=fft(hatu0.^2) #for N(u)
hatu[:,1]=fft(hatu0) #first col of hatu is fft of IC
hatNu=G.*fft(hatu2) #N(u)=u u_x = (1/2)D_x(u^2) in Fourier basis
dt2=dt/2
dt3=3*dt/2
B=(ones(Nx)-dt2*L).^-1
A=ones(Nx)+dt2*L
for n=1:Nt-1
hatNu1=hatNu
hatNu=G.*(fft(real(ifft(hatu[:,n]))).^2)
hatu[:,n+1]=B.*(A.*hatu[:,n]+dt3*hatNu-dt2*hatNu1)
usol[:,n]= real(ifft(hatu[:,n]))#
end
end
Out[4]:
In [9]:
#CNAB time-stepping for Kuramoto-Sivashinsky Equation
# u_t(x,t)=-u_{xx}(x,t)-u_{xxxx}(x,t)-(u u_x)(x,t)
# x in [0,Lx], discretized into Nx intervals of size Lx/Nx
# t in [0,T], discretized into T/dt intervals of size dt
using PyPlot
#parameter definitions
Lx=100
t=100
Nx=256
dt=1/16
dx=float(Lx/Nx)
x=collect(0:dx:Lx-dx)
T=collect(0:dt:t-dt)
Nt=round(Int,t./dt)
a=2.5*pi/Lx
u=cos(a*x)+.01*cos(2*a*x)
function cnabks(u,Lx,T,Nx,dt)
kx=vcat(0:Nx/2-1,0,-Nx/2+1:-1) #integer wavenumbers
alpha=2*pi*kx/Lx #non-dimensionalized wavenumbers
#operators in Fourier basis
D=1im*alpha #D_x(e^{i alpha x}=i alpha e^{i alpha x}
L=-D.^2-D.^4 #Lu=-u_{xx}-u_{xxxx}
G=-.5*D #(1/2)D_x u^2= u u_x
hatu=zeros(Complex{Float64},Nx,Nt)
hatNu=zeros(Complex{Float64},Nx,Nt)
usol=zeros(Complex{Float64},Nx,Nt)
hatu0=fft(u)
hatu2=fft(hatu0.^2) #for N(u)
hatu[:,1]=fft(hatu0) #first col of hatu is fft of IC
hatNu=G.*fft(hatu2) #N(u)=u u_x = (1/2)D_x(u^2) in Fourier basis
dt2=dt/2
dt3=3*dt/2
B=(ones(Nx)-dt2*L).^-1
A=ones(Nx)+dt2*L
for n=1:Nt-1
hatNu1=hatNu
hatNu=G.*(fft(real(ifft(hatu[:,n]))).^2)
hatu[:,n+1]=B.*(A.*hatu[:,n]+dt3*hatNu-dt2*hatNu1)
usol[:,n]= real(ifft(hatu[:,n]))
end
return usol
end
Out[9]:
In [10]:
cnabks(u,Lx,T,Nx,dt)
Out[10]:
Revise your code from above to simulate the Swift-Hohenberg equation $u_t = (r-1) u - 2 u_{xx} - u_{xxxx} - u^3$ on a periodic spatial domain $x \in [0,L]$ with periodic boundary conditions, at parameter value $r=0.2$. Use the same discretization methods, parameters, and initial condition as in problem 1, and make the same plot.
In [7]:
#CNAB time-stepping for 1D Cubic Swift-Hohenberg Equation
# u_t(x,t)=-u_{xx}(x,t)-u_{xxxx}(x,t)-(u u_x)(x,t)
# x in [0,Lx], discretized into Nx intervals of size Lx/Nx
# t in [0,T], discretized into T/dt intervals of size dt
using PyPlot
#parameter definitions
Lx=100
t=100
Nx=256
dt=1/16
dx=float(Lx/Nx)
x=collect(0:dx:Lx-dx)
T=collect(0:dt:t-dt)
Nt=round(Int,t./dt)
a=2.5*pi/Lx
u=cos(a*x)+.01*cos(2*a*x)
function cnabsh(u,Lx,T,Nx,dt)
kx=vcat(0:Nx/2-1,0,-Nx/2+1:-1) #integer wavenumbers
alpha=2*pi*kx/Lx #non-dimensionalized wavenumbers
#operators in Fourier basis
D=1im*alpha #D_x(e^{i alpha x}=i alpha e^{i alpha x}
L=(r-1)-D.^2-D.^4 #Lu=(r-1)-u_{xx}-u_{xxxx}
hatu=zeros(Complex{Float64},Nx,Nt)
hatNu=zeros(Complex{Float64},Nx,Nt)
usol=zeros(Complex{Float64},Nx,Nt)
hatu0=fft(u)
hatu[:,1]=fft(hatu0) #first col of hatu is fft of IC
hatNu=G.*fft(hatu2) #N(u)=u u_x = (1/2)D_x(u^2) in Fourier basis
dt2=dt/2
dt3=3*dt/2
B=(ones(Nx)-dt2*L).^-1
A=ones(Nx)+dt2*L
for n=1:Nt-1
hatNu1=hatNu
hatNu=(fft(real(ifft(hatu[:,n]))).^3)
hatu[:,n+1]=B.*(A.*hatu[:,n]+dt3*hatNu-dt2*hatNu1)
usol[:,n]= real(ifft(hatu[:,n]))#
end
pcolormesh(usol)
end
Out[7]:
In [ ]: