# -*- coding: utf-8 -*-
"""
Created on Wed Apr 1 15:41:05 2020
@author: Mr.Yang
"""
"""
函数形如
Uxx+Uyy=f(x) 狄利克雷条件 步长不一样的时候 h1 != h2 任意的时候
▽•△U=f(x) 注意符号 左边是没有 - 的
U(a)=A
U(b)=B
"""
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import scipy.linalg as la
import time
#---------------------------------------------------------------------
#定义函数 先定义后调用
def F_f(x,y):
#f=0
#f=-(np.pi*np.pi-1)*np.exp(x)*np.sin(np.pi*y) #方程 1
#f=-2*np.pi*np.pi*np.sin(np.pi*x)*np.sin(np.pi*y) #方程2
#f=1.25*np.exp(x+y/2) #作业
f=-1 #测试有限元
return(f)
def e_fn(x,y):
#f1=np.exp(x)*np.sin(np.pi*y) #方程1
#f1=np.sin(np.pi*x)*np.sin(np.pi*y) #方程2
#f1=np.exp(x+y/2) #作业
f=0 #测试有限元
return(f)
#-------------------------------------------------------------------------
#-------------------------------------------------------------------------
def U_f(x,y): #上边界
#f=0 #方程 1
#f=0 #方程2
#f=np.exp(x+1/2) #作业
f=0 #测试有限元
return(f)
def D_f(x,y): #下边界
#f=0 #方程 1
#f=0 # 方程2
#f=np.exp(x+0) #作业
f=0 #测试有限元
return(f)
def L_f(x,y): #左边界
#f=np.sin(np.pi*y) #方程 1
#f=0 #方程2
#f=np.exp(0+y/2) #作业
f=0 #测试有限元
return(f)
def R_f(x,y): #右边界
#f=np.exp(1)*np.sin(np.pi*y) #方程 1
#f=0 #方程2
#f=np.exp(1+y/2) #作业
f=0 #测试有限元
return(f)
#--------------------------------------------------------------------------------
def Fdm_cd_1d(a,b,n):
m_a=a
m_b=b
m=n #区间【a,b】之间有10个点即有U1 U2 U3 U4 U5 ...U8 U9 U10
h=1/(m+1)
x = np.linspace ( m_a, m_b, m +2) #生成等差数列 x=jh
#第一项为 m_a
#最后一项为 m_b
#共有 m+2 项
#m1=100
m1=m
h1=1/(m1+1)
y = np.linspace ( m_a, m_b, m1 +2)
#print(x)
#系数矩阵A的赋值处理方法如下:
T=np.zeros((m,m))
b=np.zeros((m*m1))
#print(b)
#print('原始矩阵T如下\n')
#print(T)
for l in range(0,m-1):
#for i in range(0,2):
T[l][l] = T[l][l] -2
T[l][l+1]=T[l][l+1]+1
T[l+1][l]=T[l+1][l]+1
T[l+1][l+1]=T[l+1][l+1]+0
print('赋值后的矩阵T如下\n')
T[m-1][m-1]=-2
print(T)
F1=np.zeros((m*m1, m*m1))
for j in range(0,m1-1):
F1[j*m:(j+1)*m,j*m:(j+1)*m]=F1[j*m:(j+1)*m,j*m:(j+1)*m]+T
F1[j*m:(j+1)*m,(j+1)*m:(j+2)*m]=F1[j*m:(j+1)*m,(j+1)*m:(j+2)*m]+np.zeros((m,m))
F1[(j+1)*m:(j+2)*m,j*m:(j+1)*m]=F1[(j+1)*m:(j+2)*m,j*m:(j+1)*m]+np.zeros((m,m))
F1[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]=F1[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]+np.zeros((m,m))
F1[m*m1-m:m*m1,m*m1-m:m*m1]=T
#print(A2)
print(F1) # 得到 左边 1/(h**2)的矩阵 还有右边
F2=np.zeros((m*m1, m*m1))
for j in range(0,m1-1):
F2[j*m:(j+1)*m,j*m:(j+1)*m]=F2[j*m:(j+1)*m,j*m:(j+1)*m]+np.eye(m,m)*(-2)
F2[j*m:(j+1)*m,(j+1)*m:(j+2)*m]=F2[j*m:(j+1)*m,(j+1)*m:(j+2)*m]+np.eye(m,m)
F2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]=F2[(j+1)*m:(j+2)*m,j*m:(j+1)*m]+np.eye(m,m)
F2[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]=F2[(j+1)*m:(j+2)*m,(j+1)*m:(j+2)*m]+np.zeros((m,m))
F2[m*m1-m:m*m1,m*m1-m:m*m1]=np.eye(m,m)*(-2)
#print(A2)
print(F2)
print(F1+F2)
A2=F1/(h**2)+F2/(h**2)
print('最终所得的矩阵A如下:\n')
print(A2)
#方程右端系数矩阵b的赋值处理方法如下
# make a mistake:can't multiply sequence by non int of type float
#in R_f(x[m+1],[m]) no y
for k in range(0,m1):
for j in range(k*m,(k+1)*m):
if(k==0): #第一层的 b
if(j==k*m):
b[j]=F_f(x[j-(k*m-1)],y[k+1])-L_f(x[0],y[k+1])/(h**2)-D_f(x[1],y[0])/(h1**2)
elif(j==m*k+m-1):
b[j]=F_f(x[j-(k*m-1)],y[k+1])-R_f(x[m+1],y[k+1])/(h**2)-D_f(x[m],y[0])/(h1**2)
else:
b[j]=F_f(x[j-(k*m-1)],y[k+1])-D_f(x[j-m*k],y[0])/(h**2) #这个地方有问题 上下边值是否一样 不一样需要修改
elif(k==m-1): #最后一层的b
if(j==k*m):
b[j]=F_f(x[j-(k*m-1)],y[k+1])-L_f(x[0],y[k+1])/(h**2)-U_f(x[1],y[m+1])/(h1**2)
elif(j==m*k+m-1):
b[j]=F_f(x[j-(k*m-1)],y[k+1])-R_f(x[m+1],y[k+1])/(h**2)-U_f(x[m],y[m+1])/(h1**2)
else:
b[j]=F_f(x[j-(k*m-1)],y[k+1])-U_f(x[j-m*k],y[m+1])/(h1**2) #4月5号 由 D_f -.>U_F()
else: #中间几行的b
if(j==k*m):
b[j]=F_f(x[j-(k*m-1)],y[k+1])-L_f(x[0],y[k+1])/(h**2)
elif(j==k*m+m-1):
b[j]=F_f(x[j-(k*m-1)],y[k+1])-R_f(x[m+1],y[k+1])/(h**2)
else:
b[j]=F_f(x[j-(k*m-1)],y[k+1])
#print('最初的矩阵b如下')
#print(b)
#b=b*h*h
print('最终所得的矩阵b如下:\n')
print(b)
b=b.T
u = la.solve ( A2, b)
#print('最初的矩阵b如下')
#print(b)
#b=b*h*h
print('最终所得的矩阵b如下:\n')
print(b)
b=b.T
u = la.solve ( A2, b)
#b1=b.T 转置与否,结果 都一样
print('------------------------------------------------------------')
print('数值解计算结果u顺序依次如下:')
print(u)
print('------------------------------------------------------------')
print('数值结果排列如下:')
print ( '' )
print('U11 U12 U13 U14...')
print('U21 U22 U23 U24...')
print('U31 U32 U33 U34...')
U=u.reshape([m1,m,])
U=U.T # 需要转置 求解结果顺序为 U11 U21 U31...,U12 U22 U32...,U13 U23 U33..
print('------------------------------------------------------------')
print('即所计算的(%dx%d)个数值解:'%(m,m1))
print(U)
'''
#------------------------------------------------------------------------------------------
#精确解的3D图形程序如下
fig=plt.figure()
ax=Axes3D(fig)
x1,y1=np.meshgrid(x,y)
ax.plot_surface(x1,y1,e_fn(x1,y1),rstride=1,cstride=1,cmap='rainbow',alpha=0.9) #绘面
ax.set_xlabel('x')
ax.set_ylabel('y')
ax.set_zlabel('Real')
ax.set_title('Figure Real')
ax.grid(False)
#--------------------------------------------------------------------------
#近似解的3D图形程序如下
print ( '' )
print('------------------------------------------------------------')
print('精确解结果如下:')
print(e_fn(x1,y1).T)
# 为了画出 数值解的 三维形状图 将程序中所求出的 数值mXm的结果替换 精确解中 中心位置mxm的结果 因为边值条件是一样的
P_U= np.zeros((m+2,m1+2)) # 定义完整的数值解
#P_U1= np.zeros((m+2,m+2)) #完整的数值解
P_U=e_fn(x1,y1).T #先把精确解 赋值 给他
'''
#could not broadcast input array from shape (3,5) into shape (3,4)
'''
print(P_U[1:m+1,1:m1+1])
P_U[1:m+1,1:m1+1]=U #采取 内部 用数值解 替换 精确解 得到(m+2)X(m+2)的所有点的数值结果
print ( '' )
print('------------------------------------------------------------')
print('完整的数值解如下:')
print(P_U)
#ax.view_init(elev=300,azim=300) #旋转曲面 参数1: 从哪个高度看曲线 参数2:旋转角度
fig1=plt.figure()
ax1=Axes3D(fig1)
#x2,y2=np.meshgrid(x,y)
Z1=np.array(list(P_U),dtype='float64')
ax1.plot_surface(x1,y1,Z1.T,rstride=1,cstride=1,cmap='rainbow',alpha=0.9) #绘面
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('Numerical')
ax1.set_title('Figure Numerical')
ax1.grid(False)
#------------------------------------------------------------------------------
#误差分析 结果以及 三维图形 程序如下
err=np.zeros((m+2,m1+2))
err=abs(e_fn(x1,y1).T-P_U)
print ( '' )
print('------------------------------------------------------------')
print('误差结果如下:')
print(err)
fig1=plt.figure()
ax2=Axes3D(fig1)
#x2,y2=np.meshgrid(x,y)
Z2=np.array(list(err),dtype='float64')
surf=ax2.plot_surface(x1,y1,Z2.T,rstride=1,cstride=1,cmap='rainbow',alpha=0.9) #绘面
ax2.set_xlabel('x')
ax2.set_ylabel('y')
ax2.set_zlabel('Error')
ax2.set_title('Figure error')
#ax2.set_zlim(0,0.1)
ax2.grid(False)
cbar=plt.colorbar(surf,shrink=0.5,aspect=5)
#cbar.set_ticks([.01,.02,.03,.04,.05,.06,.07,.08,.09,.1])
#------------------------------------------------------------------------------------------
# Compare the solution and the error at the some nodes.
#
print('------------------------------------------------------------')
print ( '' )
print('只打印第三行的结果,步长为:h=%g' %(1/(m+1)))
print ( ' Node Num True Error' )
print ( '' )
t=e_fn(x1,y1).T
for i in range ( 0, m1 + 2 ):
print ( ' (%d,%d) %14.6g %14.6g %14.6g' % ( 2,i, P_U[2][i], t[2][i], err[2][i] ) )
'''
a,b,n= map(int,input('请输入左端点a,右端点b,单元格个数m:等3个系数:各数字之间空格键隔开\n').split())
start = time.clock() #输入系数 Enter 后开始计时
Fdm_cd_1d(a,b,n)
end = time.clock()
print('计算程序运行的时间\n')
print(end-start)
print('')
2D-Possion方程的FDM-齐次边界-2
最新推荐文章于 2024-09-04 10:31:20 发布
![](https://img-home.csdnimg.cn/images/20240711042549.png)