針對neumann邊界條件的差分法程式碼
需要新增人工的Dirichlet邊界
考慮的問題仍然是矩形區域的泊松方程。
如果整個區域都採用Neumann邊界條件，即
∂
Ω
N
=
∂
Ω
\partial\Omega_{N}=\partial\Omega
∂ΩN=∂Ω，則這個問題沒有定解，因為對於某個精確解
u
u
u而言，針對任意的常數C，
u
+
c
u+c
u+c都是解
為此，我們常常需要新增一個人工提供的Dirichlet邊界條件，比如說在這裡，我們選取
{
−
1
}
×
[
−
1
,
1
]
\{-1\}\times[-1,1]
{−1}×[−1,1]為Dirichlet邊界，也就是說矩形區域的左邊為Dirichlet邊界條件。那麼問題就修改為：
假設我們已經對區域進行了對應的網格剖分：
x
i
=
i
∗
δ
x
,
i
=
0
,
…
,
M
−
1
x_{i}=i*\delta x,i=0,\ldots,M-1
xi=i∗δx,i=0,…,M−1
y
j
=
j
∗
δ
y
,
i
=
0
,
…
,
N
−
1
y_{j}=j*\delta y,i=0,\ldots,N-1
yj=j∗δy,i=0,…,N−1
下面開始重點闡述在邊界上偏導數的處理。
1：一階差分近似
u N − 1 , j − u N − 2 , j δ x = u x ( x N − 1 , y j ) \frac{u_{N-1,j} - u_{N-2,j}}{\delta x}=u_{x}(x_{N-1},y_{j}) δxuN−1,j−uN−2,j=ux(xN−1,yj)，精度是 O ( h ) O(h) O(h)，其他兩條Neumann邊界處理一樣。
# -*- coding: utf-8 -*-
"""
Created on Sat Oct 10 14:09:35 2020
@author: 2001213226
"""
import torch
import time
import numpy as np
import matplotlib.pyplot as plt
import torch.nn as nn
def UU(X, order,prob):
if prob==1:
temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
if order[0]==0 and order[1]==0:
return torch.log(temp)
if order[0]==1 and order[1]==0:
return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
if order[0]==0 and order[1]==1:
return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
if order[0]==2 and order[1]==0:
return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
+ temp**(-1) * (22)
if order[0]==1 and order[1]==1:
return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
* (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
+ temp**(-1) * (18)
if order[0]==0 and order[1]==2:
return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
+ temp**(-1) * (22)
if prob==2:
if order[0]==0 and order[1]==0:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==1 and order[1]==0:
return (3*X[:,0]*X[:,0]-1) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==0 and order[1]==1:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
(torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
if order[0]==2 and order[1]==0:
return (6*X[:,0]) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==1 and order[1]==1:
return (3*X[:,0]*X[:,0]-1) * \
(torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
if order[0]==0 and order[1]==2:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if prob==3:
temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
if order[0]==0 and order[1]==0:
return temp1 * temp2**(-1)
if order[0]==1 and order[1]==0:
return (2*X[:,0]) * temp2**(-1) + \
temp1 * (-1)*temp2**(-2) * (2*X[:,0])
if order[0]==0 and order[1]==1:
return (-2*X[:,1]) * temp2**(-1) + \
temp1 * (-1)*temp2**(-2) * (2*X[:,1])
if order[0]==2 and order[1]==0:
return (2) * temp2**(-1) + \
2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
temp1 * (-1)*temp2**(-2) * (2)
if order[0]==1 and order[1]==1:
return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
(-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
if order[0]==0 and order[1]==2:
return (-2) * temp2**(-1) + \
2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
temp1 * (-1)*temp2**(-2) * (2)
def f(prob,X):
return -UU(X,[2,0],prob)- UU(X,[0,2],prob)
class FD():
def __init__(self,bound,hx,prob):
self.prob = prob
self.dim = 2
self.hx = hx
self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
self.size = self.nx[0]*self.nx[1]
self.X = torch.zeros(self.size,self.dim)
m = 0
for i in range(self.nx[0]):
for j in range(self.nx[1]):
self.X[m,0] = bound[0,0] + i*self.hx[0]
self.X[m,1] = bound[1,0] + j*self.hx[1]
m = m + 1
self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
def matrix(self):
self.A = torch.zeros(self.nx[0]*self.nx[1],self.nx[0]*self.nx[1])
dx = self.hx[0];dy = self.hx[1]
for i in range(self.nx[0]):
for j in range(self.nx[1]):
dx = self.hx[0];dy = self.hx[1]
if i == 0:
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
elif i == self.nx[0] - 1:
self.A[i*self.nx[1]+j,(i - 1)*self.nx[1]+j] = -1
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
elif j == 0:
self.A[i*self.nx[1]+j,i*self.nx[1]+j + 1] = 1
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = -1
elif j == self.nx[1] - 1:
self.A[i*self.nx[1]+j,i*self.nx[1]+j - 1] = -1
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
return self.A
def right(self):
self.b = torch.zeros(self.nx[0]*self.nx[1],1)
for i in range(self.nx[0]):
for j in range(self.nx[1]):
dx = self.hx[0];dy = self.hx[1]
X = self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:]
if i == 0:
self.b[i*self.nx[1]+j] = UU(X,[0,0],self.prob)
elif i == self.nx[0] - 1:
self.b[i*self.nx[1]+j] = UU(X,[1,0],self.prob)*dx
elif j == 0:
self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*dy
elif j == self.nx[1] - 1:
self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*dy
elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy
return self.b
def solve(self):
A = self.matrix()
b = self.right()
u,lu = torch.solve(b,A)
return u
def error(u_pred,u_acc):
temp = ((u_pred - u_acc)**2).sum()/(u_acc**2).sum()
return temp**(0.5)
bound = torch.tensor([[0,2],[0,1]]).float()
hx = [0.2,0.1]
prob = 2
fd = FD(bound,hx,prob)
u_pred = fd.solve()
u_acc = fd.u_acc
print(error(u_pred,u_acc))
這個差分法的效果極為不佳。
二階差分近似
這裡需要注意在
u
y
(
x
,
y
0
)
,
u
x
(
x
,
y
N
−
1
)
u_{y}(x,y_{0}),u_{x}(x,y_{N-1})
uy(x,y0),ux(x,yN−1)這兩個地方做二階近似的時候，係數相差一個正負號，自己最好推導一遍。
# -*- coding: utf-8 -*-
"""
Spyder Editor
This is a temporary script file.
"""
import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn
def UU(X, order,prob):
if prob==1:
temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
if order[0]==0 and order[1]==0:
return torch.log(temp)
if order[0]==1 and order[1]==0:
return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
if order[0]==0 and order[1]==1:
return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
if order[0]==2 and order[1]==0:
return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
+ temp**(-1) * (22)
if order[0]==1 and order[1]==1:
return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
* (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
+ temp**(-1) * (18)
if order[0]==0 and order[1]==2:
return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
+ temp**(-1) * (22)
if prob==2:
if order[0]==0 and order[1]==0:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==1 and order[1]==0:
return (3*X[:,0]*X[:,0]-1) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==0 and order[1]==1:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
(torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
if order[0]==2 and order[1]==0:
return (6*X[:,0]) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==1 and order[1]==1:
return (3*X[:,0]*X[:,0]-1) * \
(torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
if order[0]==0 and order[1]==2:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if prob==3:
temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
if order[0]==0 and order[1]==0:
return temp1 * temp2**(-1)
if order[0]==1 and order[1]==0:
return (2*X[:,0]) * temp2**(-1) + \
temp1 * (-1)*temp2**(-2) * (2*X[:,0])
if order[0]==0 and order[1]==1:
return (-2*X[:,1]) * temp2**(-1) + \
temp1 * (-1)*temp2**(-2) * (2*X[:,1])
if order[0]==2 and order[1]==0:
return (2) * temp2**(-1) + \
2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
temp1 * (-1)*temp2**(-2) * (2)
if order[0]==1 and order[1]==1:
return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
(-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
if order[0]==0 and order[1]==2:
return (-2) * temp2**(-1) + \
2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
temp1 * (-1)*temp2**(-2) * (2)
def f(prob,X):
return -UU(X,[2,0],prob)- UU(X,[0,2],prob)
class FD():
def __init__(self,bound,hx,prob):
self.prob = prob
self.dim = 2
self.hx = hx
self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
self.size = self.nx[0]*self.nx[1]
self.X = torch.zeros(self.size,self.dim)
m = 0
for i in range(self.nx[0]):
for j in range(self.nx[1]):
self.X[m,0] = bound[0,0] + i*self.hx[0]
self.X[m,1] = bound[1,0] + j*self.hx[1]
m = m + 1
self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
def matrix(self):
self.A = torch.zeros(self.size,self.size)
dx = self.hx[0];dy = self.hx[1]
for i in range(self.nx[0]):
for j in range(self.nx[1]):
if i == 0:
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
elif i == self.nx[0] - 1:
self.A[i*self.nx[1]+j,(i - 2)*self.nx[1]+j] = 1
self.A[i*self.nx[1]+j,(i - 1)*self.nx[1]+j] = -4
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 3
elif j == 0:
self.A[i*self.nx[1]+j,i*self.nx[1]+j + 2] = -1
self.A[i*self.nx[1]+j,i*self.nx[1]+j + 1] = 4
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = -3
elif j == self.nx[1] - 1:
self.A[i*self.nx[1]+j,i*self.nx[1]+j - 2] = 1
self.A[i*self.nx[1]+j,i*self.nx[1]+j - 1] = -4
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 3
elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
return self.A
def right(self):
self.b = torch.zeros(self.size,1)
dx = self.hx[0];dy = self.hx[1]
for i in range(self.nx[0]):
for j in range(self.nx[1]):
X = self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:]
if i == 0:
self.b[i*self.nx[1]+j] = UU(X,[0,0],self.prob)
#self.b[i*self.nx[1]+j] = UU(X,[1,0],self.prob)*2*dx
elif i == self.nx[0] - 1:
self.b[i*self.nx[1]+j] = UU(X,[1,0],self.prob)*2*dx
elif j == 0:
self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*2*dy
elif j == self.nx[1] - 1:
self.b[i*self.nx[1]+j] = UU(X,[0,1],self.prob)*2*dy
elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy
return self.b
def solve(self):
A = self.matrix()
b = self.right()
u,lu = torch.solve(b,A)
return u
def error(u_pred,u_acc):
temp = ((u_pred - u_acc)**2).sum()/(u_acc**2).sum()
return temp**(0.5)
bound = torch.tensor([[0,2],[0,1]]).float()
hx = [0.2,0.1]
prob = 3
fd = FD(bound,hx,prob)
u_pred = fd.solve()
u_acc = fd.u_acc
print(error(u_pred,u_acc))
3：構造虛擬點
這種近似方法的誤差也是
O
(
h
2
)
O(h^{2})
O(h2)，但是構造特別麻煩，在矩形區域的4個頂點上需要和其他Neumann邊界進行重新劃分
import numpy as np
import matplotlib.pyplot as plt
import torch
import time
import torch.nn as nn
def UU(X, order,prob):
if prob==1:
temp = 10*(X[:,0]+X[:,1])**2 + (X[:,0]-X[:,1])**2 + 0.5
if order[0]==0 and order[1]==0:
return torch.log(temp)
if order[0]==1 and order[1]==0:
return temp**(-1) * (20*(X[:,0]+X[:,1]) + 2*(X[:,0]-X[:,1]))
if order[0]==0 and order[1]==1:
return temp**(-1) * (20*(X[:,0]+X[:,1]) - 2*(X[:,0]-X[:,1]))
if order[0]==2 and order[1]==0:
return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) ** 2 \
+ temp**(-1) * (22)
if order[0]==1 and order[1]==1:
return - temp**(-2) * (20*(X[:,0]+X[:,1])+2*(X[:,0]-X[:,1])) \
* (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) \
+ temp**(-1) * (18)
if order[0]==0 and order[1]==2:
return - temp**(-2) * (20*(X[:,0]+X[:,1])-2*(X[:,0]-X[:,1])) ** 2 \
+ temp**(-1) * (22)
if prob==2:
if order[0]==0 and order[1]==0:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==1 and order[1]==0:
return (3*X[:,0]*X[:,0]-1) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==0 and order[1]==1:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
(torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
if order[0]==2 and order[1]==0:
return (6*X[:,0]) * \
0.5*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if order[0]==1 and order[1]==1:
return (3*X[:,0]*X[:,0]-1) * \
(torch.exp(2*X[:,1])-torch.exp(-2*X[:,1]))
if order[0]==0 and order[1]==2:
return (X[:,0]*X[:,0]*X[:,0]-X[:,0]) * \
2*(torch.exp(2*X[:,1])+torch.exp(-2*X[:,1]))
if prob==3:
temp1 = X[:,0]*X[:,0] - X[:,1]*X[:,1]
temp2 = X[:,0]*X[:,0] + X[:,1]*X[:,1] + 0.1
if order[0]==0 and order[1]==0:
return temp1 * temp2**(-1)
if order[0]==1 and order[1]==0:
return (2*X[:,0]) * temp2**(-1) + \
temp1 * (-1)*temp2**(-2) * (2*X[:,0])
if order[0]==0 and order[1]==1:
return (-2*X[:,1]) * temp2**(-1) + \
temp1 * (-1)*temp2**(-2) * (2*X[:,1])
if order[0]==2 and order[1]==0:
return (2) * temp2**(-1) + \
2 * (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,0])**2 + \
temp1 * (-1)*temp2**(-2) * (2)
if order[0]==1 and order[1]==1:
return (2*X[:,0]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
(-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,0]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,0]) * (2*X[:,1])
if order[0]==0 and order[1]==2:
return (-2) * temp2**(-1) + \
2 * (-2*X[:,1]) * (-1)*temp2**(-2) * (2*X[:,1]) + \
temp1 * (2)*temp2**(-3) * (2*X[:,1])**2 + \
temp1 * (-1)*temp2**(-2) * (2)
def f(prob,X):
return -UU(X,[2,0],prob)- UU(X,[0,2],prob)
class FD():
def __init__(self,bound,hx,prob):
self.prob = prob
self.dim = 2
self.hx = hx
self.nx = [int((bound[0,1] - bound[0,0])/self.hx[0]) + 1,int((bound[1,1] - bound[1,0])/self.hx[1]) + 1]
self.size = self.nx[0]*self.nx[1]
self.X = torch.zeros(self.size,self.dim)
m = 0
for i in range(self.nx[0]):
for j in range(self.nx[1]):
self.X[m,0] = bound[0,0] + i*self.hx[0]
self.X[m,1] = bound[1,0] + j*self.hx[1]
m = m + 1
self.u_acc = UU(self.X,[0,0],prob).view(-1,1)
def matrix(self):
self.A = torch.zeros(self.size,self.size)
dx = self.hx[0];dy = self.hx[1]
for i in range(self.nx[0]):
for j in range(self.nx[1]):
if i == 0:
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 1
elif (i == self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
elif (j == 0 and i > 0 and i < self.nx[0] - 1):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -2*dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
elif (j == self.nx[1] - 1 and i > 0 and i < self.nx[0] - 1):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -2*dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
elif (i == self.nx[0] - 1 and j == 0):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -2*dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
elif (i == self.nx[0] - 1 and j == self.nx[1] - 1):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -2*dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -2*dy/dx
elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.A[i*self.nx[1]+j,i*self.nx[1]+j] = 2*(dx/dy + dy/dx)
self.A[i*self.nx[1]+j,i*self.nx[1]+j-1] = -dx/dy
self.A[i*self.nx[1]+j,i*self.nx[1]+j+1] = -dx/dy
self.A[i*self.nx[1]+j,(i-1)*self.nx[1]+j] = -dy/dx
self.A[i*self.nx[1]+j,(i+1)*self.nx[1]+j] = -dy/dx
return self.A
def right(self):
self.b = torch.zeros(self.size,1)
dx = self.hx[0];dy = self.hx[1]
for i in range(self.nx[0]):
for j in range(self.nx[1]):
X = self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:]
if i == 0:
self.b[i*self.nx[1]+j] = UU(X,[0,0],self.prob)
elif (i == self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy
elif (j == 0 and i > 0 and i < self.nx[0] - 1):
self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy - UU(X,[0,1],self.prob)*2*dx
elif (j == self.nx[1] - 1 and i > 0 and i < self.nx[0] - 1):
self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[0,1],self.prob)*2*dx
elif (i == self.nx[0] - 1 and j == 0):
self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy - UU(X,[0,1],self.prob)*2*dx
elif (i == self.nx[0] - 1 and j == self.nx[1] - 1):
self.b[i*self.nx[1]+j] = f(self.prob,X)*dx*dy + UU(X,[1,0],self.prob)*2*dy + UU(X,[0,1],self.prob)*2*dx
elif (i > 0 and i < self.nx[0] - 1 and j > 0 and j < self.nx[1] - 1):
self.b[i*self.nx[1]+j] = f(self.prob,self.X[i*self.nx[1]+j:i*self.nx[1]+j+1,:])*dx*dy
return self.b
def solve(self):
A = self.matrix()
b = self.right()
u,lu = torch.solve(b,A)
return u
def error(u_pred,u_acc):
temp = ((u_pred - u_acc)**2).sum()/(u_acc**2).sum()
return temp**(0.5)
bound = torch.tensor([[0,2],[0,1]]).float()
hx = [0.2,0.1]
prob = 2
fd = FD(bound,hx,prob)
u_pred = fd.solve()
u_acc = fd.u_acc
print(error(u_pred,u_acc))
綜合來看，第三種方法的效果最好。
相關文章
- MySQL 針對 like 條件的優化MySql優化
- javascript對於if條件語句程式碼的優化方式JavaScript優化
- C# 針對特定的條件進行鎖操作，不用lock，而是mutexC#Mutex
- 電磁場與電磁波：讀書筆記：恆定電流邊界條件筆記
- css條紋邊框程式碼例項CSS
- 程式碼組(2)成員條件
- 技術的邊界
- 《超越邊界》
- 二分查詢左邊界，右邊界，>=,>,<=,<
- 針對高密級單位的安全保密郵件
- mysql where條件中 字串右邊的空格會忽略MySql字串
- Python 工匠：編寫條件分支程式碼的技巧Python
- loadrunner必用函式web_reg_save_param獲取多個符合邊界值條件的使用方法函式Web
- 微服務邊界微服務
- FBI針對Tor網路的惡意程式碼分析
- Dotfuscator針對C#程式碼混淆方法總結C#
- 解決多執行緒競爭條件——臨界區執行緒
- css的邊界和補白CSS
- 程式碼優化-多型代替IF條件判斷優化多型
- Golang 針對 MySQL 資料庫表結構的差異 SQL 工具GolangMySql資料庫
- 程式返回條件的0和1
- CWE-501: Trust Boundary Violation違反信任邊界的程式碼漏洞缺陷Rust
- 各瀏覽器對 onbeforeunload 事件的支援與觸發條件實現有差異瀏覽器事件
- 關於在SQL語句中ON和WHERE中條件使用的差異SQL
- css3 動態條紋邊框程式碼例項CSSS3
- 使用SSH完成條件及分頁查詢的主要程式碼
- 條件迴圈語句組成了Python程式碼的骨架Python
- 重構遺留程式碼（3）：複雜的條件語句
- 【AutoCAD .NET】如何在無邊界Hatch上選擇邊界點？
- Scala與Java差異（二）之條件控制與迴圈Java
- 這坑貨 （迭代+迭代終止條件（由誤差控制））
- [技術討論]系統間呼叫與邊界類的差別——被混淆的兩個概念
- 對於工程師責任和責任邊界的認知工程師
- oracle 12c 針對cdb的差異0備與對pdb進行恢復Oracle
- oracle中各程式的觸發條件(希望對大家有幫助) (zt)Oracle
- 關於運營邊界的思考
- MySQL版本對varchar的定義和限制條件MySql
- 針對IE6、7、8條件註釋語句(不同版本IE顯示用不用css)CSS