homework 10

code

import numpy as np
import matplotlib.pyplot as plt

def Jacobi():
    v1=np.zeros((30,30))
    a=np.linspace(11,18,8,dtype=int)
    for i in a:
        v1[i,11]=1
        v1[i,18]=-1
    v2=v1
    for k in range(200):
        for i in range(1,29):
            for j in range(1,29):
                if i in a and j==11:
                    v2[i,j]=v1[i,j]
                elif i in a and j==18:
                    v2[i,j]=v1[i,j]
                else:
                    v2[i,j]=(v1[i-1,j]+v1[i+1,j]+v1[i,j-1]+v1[i,j+1])/4
        for i in range(1,29):
            for j in range(1,29):
                if i in a and j==11:
                    v1[i,j]=v2[i,j]
                elif i in a and j==18:
                    v1[i,j]=v2[i,j]
                else:
                    v1[i,j]=(v2[i-1,j]+v2[i+1,j]+v2[i,j-1]+v2[i,j+1])/4
    return v1

   v=Jacobi()

x=np.linspace(-1,1,30)
y=np.linspace(-1,1,30)
X, Y = np.meshgrid(x, y)

plt.contour(X,Y,v, colors = 'black')
plt.show()

u,w=np.gradient(-v)
plt.quiver(X, Y, w, u)
plt.show()

ax=plt.subplot(111,projection='3d')
ax.plot_surface(X, Y, v, rstride=1, cstride=1)
plt.show()