import numpy import scipy.linalg import matplotlib.pyplot as plt L = 10.0 # Aufgabe 10.1.2 def solve_poisson1d_exact(rho, h): N = len(rho) A = numpy.zeros((N,N)) # norming elements A[0,0] = 1.0 A[N-1,N-1] = 1.0 for i in range(2,N-1): A[i,i-1] = 1.0/(h*h) for i in range(1,N-1): A[i,i] = -2.0/(h*h) for i in range(1, N-2): A[i,i+1] = 1.0/(h*h) phi = scipy.linalg.solve(A, rho) return phi # Aufgabe 10.1.3 plt.subplot(121) x = numpy.linspace(0.0, L, 1000) plt.plot(x, numpy.sin(2*numpy.pi/L*x)) plt.plot(x, -1.0/((2*numpy.pi/L)**2)*numpy.sin(2*numpy.pi/L*x)) N = 10 x = numpy.linspace(0.0, L, N) rho = numpy.sin(2*numpy.pi/L*x) phi = solve_poisson1d_exact(rho, L/N) plt.plot(x, phi, 'o-') N = 50 x = numpy.linspace(0.0, L, N) rho = numpy.sin(2*numpy.pi/L*x) phi = solve_poisson1d_exact(rho, L/N) plt.plot(x, phi, 'o-') plt.legend(('rho', 'analytical', 'N=10', 'N=50')) # Aufgabe 10.3.1 def sor_step(A, b, x, omega=1.0): N = len(b) x_new = (1.0-omega)*x for j in range(N): x_new[j] += omega/A[j,j] * \ (b[j] - numpy.dot(A[j,:j], x_new[:j]) - numpy.dot(A[j,j+1:], x[j+1:])) return x_new def sor(A, b, steps, omega=1.0): N = len(b) x = numpy.zeros_like(b) for i in range(steps): x = sor_step(A, b, x, omega) return x # Aufgabe 10.3.2 def solve_poisson1d_sor(rho, h, steps): N = len(rho) A = numpy.zeros((N,N)) # norming elements A[0,0] = 1.0 A[N-1,N-1] = 1.0 for i in range(2,N-1): A[i,i-1] = 1.0/(h*h) for i in range(1,N-1): A[i,i] = -2.0/(h*h) for i in range(1, N-2): A[i,i+1] = 1.0/(h*h) phi = sor(A, rho, steps) return phi plt.subplot(122) x = numpy.linspace(0.0, L, 1000) plt.plot(x, -1.0/((2*numpy.pi/L)**2)*numpy.sin(2*numpy.pi/L*x)) N = 10 x = numpy.linspace(0.0, L, N) rho = numpy.sin(2*numpy.pi/L*x) phi = solve_poisson1d_exact(rho, L/N) plt.plot(x, phi, 'o-') x = numpy.linspace(0.0, L, N) rho = numpy.sin(2*numpy.pi/L*x) phi = solve_poisson1d_sor(rho, L/N, 5) plt.plot(x, phi, 'o-') x = numpy.linspace(0.0, L, N) rho = numpy.sin(2*numpy.pi/L*x) phi = solve_poisson1d_sor(rho, L/N, 10) plt.plot(x, phi, 'o-') x = numpy.linspace(0.0, L, N) rho = numpy.sin(2*numpy.pi/L*x) phi = solve_poisson1d_sor(rho, L/N, 20) plt.plot(x, phi, 'o-') plt.legend(('analytical', 'exact', 'steps=5', 'steps=10', 'steps=20')) # Aufgabe plt.show()