цифровая дифференциация электрического потенциала с питоном

Question

Я использовал следующий код для получения электрического поля квадрупольного потенциала.

for n in range (nx-1): 
    for m in range (ny-1): 
    for k in range (nz-1): 
    Ex[n,m,k] =-((T[n+1,m,k]-T[n,m,k]))/((x[n+1]-x[n])); 
    Ey[n,m,k] =-((T[n,m+1,k]-T[n,m,k]))/((y[m+1]-y[m])); 
    Ez[n,m,k] =-((T[n,m,k+1]-T[n,m,k]))/((z[k+1]-z[k])); 
    return Ex,Ey,Ez, T 

здесь Т потенциал в 3D получено путем решения уравнения Лапласа численно, как вы можете видеть из фигуры красного электрода (положительными) имеет неправильный вектор электрического поля направления (правую верхнюю сторону и правую вниз сторону верхнего электрода) та же ошибка также в других электродах. то есть отрицательный электрод имеет выходящий вектор электрического поля, который должен быть в противоположном направлении.

Я также использовал центральный метод разницы, но я получаю такую ​​же цифру. Не могли бы вы рассказать мне, что не так с моей дифференциацией?

Answer 1

Ваша проблема заключается в том, каким образом Matplotlib интерпретирует направления, которое само по себе берет свое начало на пути матрицы обычно индексируются, что обратное тому, что многие люди (например, я) думаю. В частности, первым индексом является индекс вертикальный индекс, а второй индекс - , горизонтальный индекс. Это было бы ясно, если бы вы использовали неквадратный массив, который бы показал, что ваши x и y были назад. Следующий пример иллюстрирует, как ваш градиент (только слегка измененный) дает правильный результат. На первом рисунке показано, что вы на самом деле замышляли, когда компоненты x и y градиентов менялись местами. Рисунок (который вы можете получить, запустив это) показывает, что градиент не ортогонален контурным линиям (как и должно быть) и иногда идет в неправильном направлении.

На втором и третьем рисунках показан правильный способ построения градиента с использованием хорошего подхода, в котором у нас есть 2D-массивы для координат, а также объекты, которые мы рисуем (quiver или contour). Я следую вашему коду при использовании meshgrid, чтобы сгенерировать мои X и Y, но вам нужно было поменять аргументы на meshgrid, чтобы получить правильные размеры. Более чистый подход состоял бы в том, чтобы использовать один и тот же стиль кодирования для генерации координат в качестве построенной вещи, в этом случае цикл while с явными индексами был бы оптимальным.

import numpy as np 
import matplotlib.pyplot as plt 

x = np.arange(-3,3,0.2) 
y = np.arange(-5,5,0.2) 
z = np.arange(-5,5,0.2) 

def grad(f): 
    gx = np.zeros_like(f) 
    gy = np.zeros_like(f) 
    gz = np.zeros_like(f) 
    for n in range(1,len(x)-1): 
     for m in range(1,len(y)-1): 
      for k in range(1,len(z)-1): 
       gx[n,m,k] = (T[n+1,m,k]-T[n-1,m,k])/(x[n+1]-x[n-1]); 
       gy[n,m,k] = (T[n,m+1,k]-T[n,m-1,k])/(y[m+1]-y[m-1]); 
       gz[n,m,k] = (T[n,m,k+1]-T[n,m,k-1])/(z[k+1]-z[k-1]); 
    return gx, gy, gz 

T = np.zeros((len(x), len(y), len(z))) 
for n in range(len(x)): 
    for m in range(len(y)): 
     for k in range(len(z)): 
      T[n,m,k] = np.sin((x[n] - y[m])/3.0) + 0.3*np.cos(y[m]) + z[k]**2 

gx,gy,gz = grad(T) 

Y, X= np.meshgrid(y,x) 

plt.figure('WRONG with x and y') 
plt.contour(y, x, T[:,:,round(len(z)/2)], 64) 
plt.colorbar() 
plt.quiver(y, x, 10*gx[:,:,round(len(z)/2)], 10*gy[:,:,round(len(z)/2)]) 
plt.xlabel("x") 
plt.ylabel("y") 
plt.axes().set_aspect('equal') 


plt.figure('with X and Y') 
plt.contour(X, Y, T[:,:,round(len(z)/2)], 64) 
plt.colorbar() 
plt.quiver(X, Y, 10*gx[:,:,round(len(z)/2)], 10*gy[:,:,round(len(z)/2)]) 
plt.xlabel("X") 
plt.ylabel("Y") 
plt.axes().set_aspect('equal') 

plt.figure('with Y and X') 
plt.contour(Y, X, T[:,:,round(len(z)/2)], 64) 
plt.colorbar() 
plt.quiver(Y, X, 10*gy[:,:,round(len(z)/2)], 10*gx[:,:,round(len(z)/2)]) 
plt.xlabel("Y") 
plt.ylabel("X") 
plt.axes().set_aspect('equal') 

plt.show() 

Наконец, я укажу еще раз, чистый способ выявить эту ошибку было бы запустить программу с nx != ny. Код не удался с сообщением об ошибке, указывающим, что размеры массива не совпадают, что привело бы вас к замене вещей (надеюсь) правильным образом.