Numerical methods for solving elliptic equations

from numpy import*
from numpy.linalg import solve
import matplotlib.pyplot as plt
x0=0#Начальная координата области решения
xn=5#Конечная координата области решения
n=100#Число точек координатной сетки
dx=(xn-x0)/(n-1)#Задание равномерной координатной сетки с шагом dx
x=[i*dx+x0 for i in arange(0,n,1)]#Задание равномерной координатной сетки с шагом dxdeff(i):#Функция правой части уравненияreturn2*sin(x[i]**2)+cos(x[i]**2)
v1=1.0#Вид ГУ на левой границе (1 - Дирихле, 2 - Неймана)
g1=0.0#Значение ГУ на левой границе
v2=2.0#'Вид ГУ на правой границе (1 - Дирихле, 2 - Неймана)
g2=-0.5#Значение ГУ на правой границе
a=zeros([n,n])#Задание матрицы коэффициентов СЛАУ размерностью n x n
b=zeros([1,n])# Задание матрицы-строки свободных членов СЛАУ размерностью 1 x n#Определение коэффициентов и свободных членов СЛАУ,# соответствующих граничным условиям и проверка корректности#значений параметров v1, v2
b[0,n-1]=g1;
if v1==1:
         a[0,0]=1elif v1==2:
         a[0,0]=-1/dx
         a[0,1]=1/dx;
else:
         print('Параметр v1 имеет неправильное значение')
b[0,n-1]=g2;
if v2==1:
         a[n-1,n-1]=1elif v2==2:
         a[n-1,n-1]=1/dx
         a[n-1,n-2]=-1/dx;
else:
         print('Параметр v2 имеет неправильное значение')
#Определение коэффициентов и свободных членов СЛАУ,# соответствующих внутренним точкам области         for i in arange(1, n-1,1):
         a[i,i]=-2/dx**2
         a[i,i+1]=1/dx**2
         a[i,i-1]=1/dx**2
         b[0,i]=f(i)
u=linalg.solve(a,b.T).T#Решение СЛАУdefviz(v1,v2):if v1==v2==1:
                  return"ГУ  Дирихле на левой и ГУ Дирихле на правой  границе "elif v1==1and v2==2:
                  return"ГУ  Дирихле на левой и ГУ Неймана на правой  границе "elif v2==1and v2==1:
                  return"ГУ  Неймана на левой и ГУ Дирихле  на правой  границе "
plt.figure()
plt.title("График функции правой части уравнения Пуассона")
y=[f(i) for i in arange(0,n,1)]
plt.plot(x,y)
plt.grid(True)
plt.xlabel('x')
plt.ylabel('f(x)')
plt.figure()
plt.title("График искомой функции уравнения Пуассона")
plt.xlabel('x')
plt.ylabel('u(x)')
plt.plot(x,u[0,:],label='%s'%viz(v1,v2))
plt.legend(loc='best')
plt.grid(True)
plt.show()

from __future__ import print_function 
from __future__ import division 
import numpy as np 
import time 
ti = time.clock() 
m = 1000 
A = np.zeros((m, m)) 
B = np.zeros((m, 1))  
A[0, 0] = 1 
A[0, 1] = 2 
B[0, 0] = 1for i in range(1, m-1): 
    A[i, i-1] = 7 
    A[i, i] = 8  
    A[i, i+1] = 9 
    B[i, 0] = 2 
A[m-1, m-2] = 3 
A[m-1, m-1] = 4 
B[m-1, 0] = 3 
print('A \n', A) 
print('B \n', B) 
x = np.linalg.solve(A, B)  # solve A*x = B for x 
print('x \n', x) 
print('NUMPY time', time.clock()-ti, 'seconds') 

from numpy import *
""" Численное решение задачи Дирихле
для уравнения конвекции-диффузии в
прямоугольнике.Метод релаксации."""defrelaxation(b, f, I1, I2, n1, n2, omega, tol = 1.e-8):
         h1 = I1 / n1
         h2 = I2 / n2
         d = 2. / h1**2 + 2. / h2**2
         y = zeros([n1+1, n2+1])
         ff = zeros([n1+1, n2+1])
         bb = zeros([n1+1, n2+1])
         for j in arange(1,n2,1):
                  for i in arange(1,n1,1):
                           ff [i,j] = f(i*h1, j*h2)
                           bb[i,j] = b(i*h1, j*h2)                           
         #максимальное число итераций - 10000for k in arange(1, 10001,1):
                  rn = 0.for j in arange(1,n2,1):
                           for i in arange(1,n1,1):                                    
                                    rr = - (y[i-1,j] - 2.*y [i, j] + y[i+1,j]) / h1**2 \
                                         - (y[i,j-1] - 2.*y [i,j] + y[i,j+1]) / h2**2 \
                                         + bb[i,j]*(y [i+1,j] - y [i-1,j]) / (2.*h1) - ff [i,j]                                    
                                    rn = rn + rr**2
                                    y[i,j] = y[i,j] - omega * rr / d
                  rn = rn*h1*h2
                  if rn < tol**2: return y, k                  
         print ('Метод релаксации не сходиться:')
         print ('после 10000 итерации остаток=',sqrt(rn))
import matplotlib.pyplot as plt
bcList = [0., 10.]
sglist = ['-','--']
kk = 0for bc in bcList:         
         I1 = 1.
         I2 = 1.deff(x,y):return1.defb(x,y):return bc
         n1 = 25
         n2 = 25
         m = 20
         om = linspace(1., 1.95, m)
         it = zeros(([m]))
         for k in arange(0,m,1):
                  omega = om[k]
                  y, iter = relaxation(b, f, I1, I2, n1, n2, omega, tol=1.e-6)
                  it[k] = iter
         s1= 'b =' + str(bc)
         sg = sglist[kk]
         kk = kk+1
         plt.plot( om,it, sg, label = s1)
plt.title("Число итераций метода релаксации\n для приближённого решения эллиптической задачи\n с использованием заданного параметра релаксации $\\omega$")
plt.xlabel('$\\omega$')
plt.ylabel('iterations')
plt.legend(loc=0)
plt.grid(True)
plt.show(