In [1]:
3+5*6
Out[1]:
33
In [2]:
import math
math.sqrt(3**2+4**2)
Out[2]:
5.0
In [3]:
1/7
Out[3]:
0.14285714285714285
In [4]:
a=3+4j
print(a)

(3+4j)
In [5]:
a.real
Out[5]:
3.0
In [6]:
math.exp(1)**math.pi-math.pi**math.exp(1)
Out[6]:
0.6815349144182221
In [7]:
abs(a)
Out[7]:
5.0
In [8]:
print([a.real, a.imag, abs(a)])

[3.0, 4.0, 5.0]
In [9]:
abs(math.cos(math.pi/6)+math.sin(math.pi/6)*1j)
Out[9]:
1.0
In [10]:
import numpy as np
a=np.array([1,2,3])
print(a)

[1 2 3]
In [11]:
np.log(a)
Out[11]:
array([ 0.        ,  0.69314718,  1.09861229])
In [12]:
np.power(a,2)
Out[12]:
array([1, 4, 9])
In [13]:
a=np.linspace(1,5,3)
print(a)

[ 1.  3.  5.]
In [14]:
a=np.linspace(0,10.0,10)
print(a)

[  0.           1.11111111   2.22222222   3.33333333   4.44444444
   5.55555556   6.66666667   7.77777778   8.88888889  10.        ]
In [15]:
a=np.logspace(-1,3,5)
print(a)

[  1.00000000e-01   1.00000000e+00   1.00000000e+01   1.00000000e+02
   1.00000000e+03]
In [16]:
a=np.array([1,2,3]); b=np.array([2,4,5])
print(a,b)

[1 2 3] [2 4 5]
In [17]:
a+b
Out[17]:
array([3, 6, 8])
In [18]:
a+1
Out[18]:
array([2, 3, 4])
In [19]:
a*b
Out[19]:
array([ 2,  8, 15])
In [20]:
a=np.array([[1],[2],[3]])
print(a)
print(a.T)

[[1]
 [2]
 [3]]
[[1 2 3]]
In [21]:
a=np.ones(20)
print(a)

[ 1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.  1.
  1.  1.]
In [22]:
b=np.zeros(5)
print(b)

[ 0.  0.  0.  0.  0.]
In [23]:
c=np.random.rand(5)
print(c)

[ 0.78264285  0.21923003  0.22951218  0.15023256  0.28353315]
In [24]:
print(c[0:3])

[ 0.78264285  0.21923003  0.22951218]
In [25]:
print(c[2:4])

[ 0.22951218  0.15023256]
In [26]:
c[4]
Out[26]:
0.28353314677174757
In [27]:
sum(c)
Out[27]:
1.665150771529968
In [28]:
a=np.array([1,2,3,4,5])
print(a)

[1 2 3 4 5]
In [29]:
np.cumsum(a)
Out[29]:
array([ 1,  3,  6, 10, 15])
In [30]:
np.diff(a)
Out[30]:
array([1, 1, 1, 1])
In [31]:
import matplotlib.pyplot as plt
x=np.array([1,2,3,5,7])
y=np.array([2,1,4,6,5])
plt.figure(figsize=(15,5))
plt.subplot(1,3,1)
plt.plot(x,y)
plt.subplot(1,3,2)
plt.plot(x,y,'-o')
plt.subplot(1,3,3)
plt.plot(x,y,'.')
plt.show()

<matplotlib.figure.Figure at 0x7fcca2ffe978>
In [32]:
x=np.linspace(1,6,1000)
y=math.exp(1)*np.log(x)/x
plt.plot(x,y)
plt.grid()
plt.show()
In [33]:
def Logi_map(n,a):
    x=np.zeros(n)
    
    x[0]=0.12345
    for k in range(n-1):
        x[k+1]=a*x[k]*(1-x[k])
        #print(k)
    
    plt.plot(x,'-o')
    plt.show()
In [34]:
Logi_map(101,3.5)
In [35]:
Logi_map(100,1)
In [36]:
Logi_map(100,2)
In [37]:
Logi_map(100,3)
In [38]:
Logi_map(200,4)

Henon_map

calculate the equation

$x_{k+1}=ax_k(1-x_k)$.

In [39]:
def Henon_map(n,a,b):
    x=np.zeros(n)
    y=np.zeros(n)
    
    x[0]=0
    y[0]=0.2
    for k in range(n-1):
        x[k+1]=1-a*x[k]**2+y[k]
        y[k+1]=b*x[k]
        #print(k)
    
    plt.plot(x,y,'.')
    plt.show()
In [40]:
Henon_map(10000,1.4,0.3)
In [41]:
Henon_map(20000,0.2,0.99)
In [42]:
Henon_map(10000,0.2,-0.9991)

Consider a test periodic function,

$u(x)=\sin(\exp(3\sin(x)))$.

The period is 2$\pi$.

  1. Discretize $x$ into $n$ grids, and define vectors $u$ and $xx$ as explained in the previous page.
  2. Use "centered" and "higher-order centered"shemes to evaluate $u'$.Name them $v_{c1}$ and $v_{c2}$.
    Smilarly, use "centered" and "higher-order centered" sheme to evaluate $u''$, and name them $w_{c1}$ and $w_{c2}$.
  3. Use FFT to evaluate $u'$ and $u''$. Name them $v_f$ and $w_f$, respectively.
  4. $u'$ and $u''$ can be analytically derived. First, let us rewrite $$u(x)=\sin(h)$$ with $$h(x)=\exp(3\sin(x))$$ then $$u'(x)=\cos(h)\cdot3h\cos x$$ $$u''(x)=3\cos(h)(3\cos^2x-\sin x)h-\sin(h)9h^2\cos^2 x$$ call these analytical solutions $v_e$ and $w_e$.
  5. Compare $v_{c1}, v_{c2}, v_f$ and $v_e$. Also, compare $w_{c1}, w_{c2}, w_f, $ and $w_e$.
  6. Plot $\mathrm{abs}(v_{c1}-v_e)$, $\mathrm{abs}(v_{c2}-v_e)$, $\mathrm{abs}(v_f -v_e)$. Also, plot $\mathrm{abs}(w_f-w_e)$.
In [43]:
import numpy as np
import matplotlib.pyplot as plt
def discritization(n):
    xmax=2*np.pi; dx=xmax/n
    jj=np.arange(n); xx=jj*dx
    
    h=np.exp(3*np.sin(xx)); u=np.sin(h)
    v=np.cos(h)*3*h*np.cos(xx)
    w=3*np.cos(h)*(3*np.cos(xx)**2-np.sin(xx))*h-np.sin(h)*9*h**2*np.cos(xx)**2
    
    plt.figure(figsize=(15,10))
    plt.subplot(3,2,1)
    plt.plot(xx,u,'-k')
    
    jp1=jj+1; jp1[n-1]=1
    jp2=jj+2; jp2[n-2]=1; jp2[n-1]=2
    jm1=jj-1; jm1[0]=n-1; 
    jm2=jj-2; jm2[0]=n-2; jm2[1]=n-1

    up1=u[jp1]; up2=u[jp2]; um1=u[jm1]; um2=u[jm2]
    
    #centered
    vc1=(up1-um1)/(2*dx)
    wc1=(up1-2*u+um1)/dx**2
    
    #centered (higher order)
    vc2=(-up2+8*up1-8*um1+um2)/(12*dx)
    wc2=(-up2+16*up1-30*u+16*um1-um2)/(12*dx**2)
    
    #fft
    fu=np.fft.fft(u)
    kk=2*np.pi/xmax*np.r_[np.arange(n/2),np.arange(-n/2,0,1)]
    vf=np.real(np.fft.ifft(1j*kk*fu))
    wf=np.real(np.fft.ifft(-kk**2*fu))
    #print(len(kk))
    #print(kk)
    plt.subplot(3,2,3)
    plt.plot(xx,vc1,xx,vc2)
    
    plt.subplot(3,2,5)
    plt.plot(xx,wc1)
    
    plt.subplot(3,2,4)
    plt.semilogy(xx,np.abs(vc1-v),xx,np.abs(vc2-v),xx,np.abs(vf-v))
    
    plt.subplot(3,2,6)
    plt.semilogy(xx,np.abs(wc1-w),xx,np.abs(wc2-w),xx,np.abs(wf-w))
    plt.show()
In [44]:
discritization(256)

FTCS (Forward Time Centered Space)

Consider a partial dirrential equation (PDE) for $u(x,t)$ $$\frac{\partial u}{\partial t}+c\frac{\partial u}{\partial x}=0$$ where $c$ is a constant. A simple way to discretize the equation above is to use "forward time" and "centered space" shemes, i.e., $$\frac{\partial u(x,t)}{\partial t}\rightarrow\frac{u_j^{k+1}-u_j^k}{\Delta t} \qquad\qquad\frac{\partial u(x,t)}{\partial x}\rightarrow\frac{u_{j+1}^k-u_{j-1}^k}{2\Delta x}$$ which lead to an explicit expression for $u_j^{k+1}$ $$u_j^{k+1}=u_j^k-\frac{c\Delta t}{2\Delta x}(u_{j+1}^k-u_{j-1}^k)$$.

In [45]:
#%matplotlib nbagg
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
#plt.ion()
#from IPython.display import display, clear_output
import matplotlib.animation as animation

def FTCS(dx,nx,dt,nt):
    u=np.r_[np.zeros(int(nx/4)), np.ones(int(nx/2)), np.zeros(int(nx/4))]
    
    jj=np.arange(nx); jp=np.r_[jj[0:nx-1]+1,0]; jm=np.r_[nx-1,jj[1:nx]-1] #jm=jj-1; jm[0]=nx-1
    #jj=np.arange(nx); jp=jj+1; jp[nx-1]=1; jm=jj-1; jm[0]=nx-1
    #print(jj, jp, jm)
    xx=dx*jj; xmax=dx*nx; c=1
    plt.plot(xx,u,'-o')
    plt.show()
    
    cf=c*dt/2/dx
    
    fig = plt.figure()
    ims = []
    #axe = fig.add_subplot(111)
    for it in range(nt):
        #clear_output(wait=True)
        u=u-cf*(u[jp]-u[jm])
        
        plt.ylim(-1,2)
        im=plt.plot(xx,u,'k-o')
        ims.append(im)
    ani = animation.ArtistAnimation(fig, ims)
    #plt.show()
    #ani.save('hoge.gif')
    return(ani)
In [46]:
FTCS(1,100,0.1,200)
Out[46]:
<matplotlib.animation.ArtistAnimation at 0x7fcca004a518>

Lax-Friedrichs

Let us try the different scheme, Lax-Friedrichs: $$\frac{\partial u(x,t)}{\partial t}\rightarrow\frac{u_j^{k+1}-(u_{j+1}^k+u_{j-1}^k)/2}{\Delta t}\qquad\qquad\frac{\partial u(x,t)}{\partial x}\rightarrow\frac{u_{j+1}^k-u_{j-1}^k}{2\Delta x}$$ Then, one again has an iteration relation $$u_j^{k+1}=\frac{u_{j+1}^k+u_{j-1}^k}{2}-\frac{c\Delta t}{2\Delta x}(u_{j+1}^k-u_{j-1}^k)$$

In [47]:
#%matplotlib nbagg
import matplotlib
import numpy as np
import matplotlib.pyplot as plt
plt.ion()
from IPython.display import display, clear_output
#import matplotlib.animation as animation

def LaxFried(dx,nx,dt,nt):
    u=np.r_[np.zeros(int(nx/4)), np.ones(int(nx/2)), np.zeros(int(nx/4))]
    
    jj=np.arange(nx); jp=np.r_[jj[0:nx-1]+1,0]; jm=np.r_[nx-1,jj[1:nx]-1] #jm=jj-1; jm[0]=nx-1
    #print(jj, jp, jm)
    xx=dx*jj; xmax=dx*nx; c=1.
    plt.plot(xx,u,'-o')
  #  plt.show()
    
    cf=c*dt/2/dx
    print(cf)
    fig = plt.figure()
    axe = fig.add_subplot(111)
    for it in range(nt):
        clear_output(wait=True)
     
        u=(u[jp]+u[jm])/2-cf*(u[jp]-u[jm])
        
        plt.ylim(-1,2)
        plt.xlim(0,xmax)
        axe.plot(xx,u,'-o')
        #fig.set_size_inches(18.5, 10.5)
        display(fig)
        axe.cla()
In [48]:
LaxFried(1,256,.1,100)

Matrix

In [49]:
import numpy as np
a=np.array([[2,1,3,5,4,5],[1,3,5,5,2,2],[3,3,5,3,2,2],[2,3,1,3,2,3]])
a[2,:]
a[:,2]
Out[49]:
array([3, 5, 5, 1])
In [50]:
xx,yy=np.meshgrid(np.linspace(1,3,3),np.linspace(2,5,4))
xx
yy
Out[50]:
array([[ 2.,  2.,  2.],
       [ 3.,  3.,  3.],
       [ 4.,  4.,  4.],
       [ 5.,  5.,  5.]])
In [51]:
def plots():
    x,y=np.meshgrid(np.linspace(-10,10,50),np.linspace(-10,10,50))
    r=np.sqrt(x**2+y**2)
    f=np.sin(r)/r
    #print(x)
    plt.subplot(2,2,1); plt.contour(x,y,f)
    
    plt.subplot(2,2,2); plt.pcolor(x,y,f)
    
    plt.show()
In [52]:
plots()
In [53]:
def vecplot():
    x,y=np.meshgrid(np.linspace(-4,4,18),np.linspace(-4,4,18))
    px=x; py=-y
    plt.quiver(x,y,px,py)
    plt.show()
In [54]:
vecplot()

Lax-Wendrof

In [55]:
#%matplotlib nbagg
import matplotlib
import numpy as np
import matplotlib.pyplot as plt

def LaxWend(dx,nx,dt,nt):
    u=np.r_[np.zeros(int(nx/4)), np.ones(int(nx/2)), np.zeros(int(nx/4))]
    
    jj=np.arange(nx); jp=np.r_[jj[0:nx-1]+1,0]; jm=np.r_[nx-1,jj[1:nx]-1] #jm=jj-1; jm[0]=nx-1
    #print(jj, jp, jm)
    xx=dx*jj; xmax=dx*nx; c=1; cc=c*dt/dx; cch=0.5*cc
    tt=dt*np.arange(nt)
    #plt.plot(xx,u,'-o')
  
    udata=np.zeros((nt,nx))
    udata[0,:]=u
    
    cf=c*dt/2/dx
    #print(cf)
    #fig = plt.figure()
    #axe = fig.add_subplot(111)
    for it in range(nt):
        #clear_output(wait=True)
        u=u-cch*(u[jp]-u[jm])+cc*cch*(u[jp]-2*u+u[jm])
        #u=u-cf*(u[jp]-u[jm])
        udata[it,:]=u
        #u=(u[jp]+u[jm])/2-cf*(u[jp]-u[jm])

    plt.pcolor(xx,tt,udata); plt.show()
    plt.plot(xx,u); plt.show()
In [56]:
%time LaxWend(1,256,0.1,1000)

CPU times: user 17.3 s, sys: 1.44 s, total: 18.7 s
Wall time: 20.2 s

Runge-Kutta+FFT

In [57]:
#%matplotlib nbagg
import matplotlib
import numpy as np
import matplotlib.pyplot as plt

def rk4fft(dx,nx,dt,nt):
    u=np.r_[np.zeros(int(nx/4)), np.ones(int(nx/2)), np.zeros(int(nx/4))]
    #u=np.sin(2*np.pi/nx*np.arange(nx)*2)
    jj=np.arange(nx)#; jp=np.r_[jj[0:nx-1]+1,0]; jm=np.r_[nx-1,jj[1:nx]-1] #jm=jj-1; jm[0]=nx-1
    #print(jj, jp, jm)
    xx=dx*jj; xmax=dx*nx; c=1
    tt=dt*np.arange(nt)
    kk=2*np.pi/xmax*np.r_[np.arange(nx/2),0,np.arange(-nx/2+1,0,1)]
    
    plt.plot(xx,u,'-o')
  
    udata=np.zeros((nt,nx))
    udata[0,:]=u
    
    #cf=c*dt/2/dx
    #print(cf)
    fig = plt.figure()
    axe = fig.add_subplot(111)
    for it in range(nt):
        #clear_output(wait=True)
        kk1=-c*dt*np.real(np.fft.ifft(1j*kk*np.fft.fft(u)))
        kk2=-c*dt*np.real(np.fft.ifft(1j*kk*np.fft.fft(u+0.5*kk1)))
        kk3=-c*dt*np.real(np.fft.ifft(1j*kk*np.fft.fft(u+0.5*kk2)))
        kk4=-c*dt*np.real(np.fft.ifft(1j*kk*np.fft.fft(u+    kk3)))
        u  =u+(kk1+2*kk2+2*kk3+kk4)/6
        
        #u=u-cch*(u[jp]-u[jm])+cc*cch*(u[jp]-2*u+u[jm])
        udata[it,:]=u
        #u=(u[jp]+u[jm])/2-cf*(u[jp]-u[jm])
        
    plt.pcolor(xx,tt,udata); plt.show()
    plt.plot(xx,u,'-o'); plt.show()
In [58]:
%time rk4fft(1,256,0.1,1024)

CPU times: user 16.8 s, sys: 984 ms, total: 17.8 s
Wall time: 18.7 s

Julia Map

In [59]:
def julia(c):
    size=400
    x=np.linspace(-1.5,1.5,size)
    y=np.linspace(-1.0,1.0,size)
    z=np.zeros((size,size),dtype=complex)
    julia=np.zeros((size,size))
    
    X,Y=np.meshgrid(x,y)
    z=X+Y*1j

    it=0
    while z.any!=False:

        z[z!=False]=z[z!=False]**2+c

        ind=np.where(np.abs(z)>2.0)
        julia[ind]=1/(10+it)
        z[ind]=False
        
        it+=1
        if it>80: break

    plt.figure(figsize=(20,10))
    plt.imshow(julia,origin='lower',cmap='terrain')
    plt.colorbar()
    plt.show()
In [60]:
c=-0.70176-0.3842j
%time julia(c)

CPU times: user 1.67 s, sys: 78.1 ms, total: 1.75 s
Wall time: 1.84 s
In [61]:
c=0.285+0.01j
%time julia(c)

CPU times: user 1.77 s, sys: 62.5 ms, total: 1.83 s
Wall time: 1.96 s
In [62]:
c=0.0000-0.8000j
%time julia(c)

CPU times: user 1.73 s, sys: 62.5 ms, total: 1.8 s
Wall time: 2.04 s
In [63]:
c=-0.4000+0.6000j
%time julia(c)

CPU times: user 2.12 s, sys: 141 ms, total: 2.27 s
Wall time: 2.75 s
In [64]:
c=-0.8+0.156j
%time julia(c)

CPU times: user 1.48 s, sys: 46.9 ms, total: 1.53 s
Wall time: 1.57 s