"""
Created on Tue Mar 10 14:24:05 2020
@author: Mr.Yang
"""
"""
函数形如
u''(x)=f(x)
u(a)=A
u(b)=B
例题1：u''(x)=sin(x)
u(0)=
u(1)=
例题2：u''(x)=x*x*exp(x)
u(0)=
u(1)=
f(x)=x*x*exp(x)
例题3：
The PDE is defined for 0 < x < 1:
# -u'' = f
# with right hand side
# f(x) = -(exact(x)'')
# and boundary conditions
# u(0) = exact(0),
# u(1) = exact(1).
#
# The exact solution is:解析解
# exact(x) = x * ( 1 - x ) * exp ( x )
# The boundary conditions are
# u(0) = 0.0 = exact(0.0),
# u(1) = 0.0 = exact(1.0).
# The right hand side is:
# f(x) = x * ( x + 3 ) * exp ( x ) = - ( exact''(x) )
"""
import matplotlib.pyplot as plt
import numpy as np
import scipy.linalg as la
import time
def f_x(x):
f=np.sin(np.pi*x)
return(f)
def e_fn(x):
f1=-1/np.pi**2*np.sin(np.pi*x)
return(f1)
def Fdm_cd_1d(a,b,n):
m_a=a
m_b=b
m=n
h=1/(m+1)
x = np.linspace ( m_a, m_b, m +2)
print(x)
A=np.zeros((m,m))
b=np.zeros((m))
print('原始矩阵A如下\n')
print(A)
for l in range(0,m-1):
A[l][l]=A[l][l]+0
A[l][l+1]=A[l][l+1]+1
A[l+1][l]=A[l+1][l]+1
A[l+1][l+1]=A[l+1][l+1]-2
print('赋值后的矩阵A如下\n')
print(A)
A[0][0]=-2
A=A/(h*h)
print('最终所得的矩阵A如下:\n')
print(A)
for i in range(0,m):
b[i]=b[i]+f_x(x[i+1])
print(b)
b[0]=b[0]-e_fn(m_a)/(h*h)
b[m-1]=b[m-1]-e_fn(m_b)/(h*h)
print('最终所得的矩阵b如下：\n')
print(b)
b=b.T
u = la.solve ( A, b)
print('数值解计算结果u如下：')
print(u)
npp = 51
xp = np.linspace (m_a,m_b,npp)
up = np.zeros ( npp )
for i in range ( 0, npp):
up[i] = e_fn ( xp[i] )
U=np.zeros((m+2))
for j in range(0,m+2):
if(j==0):
U[j]=e_fn(m_a)
elif(j==m+1):
U[j]=e_fn(m_b)
else:
U[j]=u[j-1]
print('')
print('已知有m个U的值，再加上左端点值，右端点值')
print(U)
uex = np.zeros ( m + 2 )
for i in range ( 0, m + 2 ):
uex[i] = e_fn(x[i] )
print ( '' )
print ( ' Node 数值解 解析解 Error' )
print ( '' )
err=np.zeros(m+2)
for i in range ( 0, m + 2 ):
err[i] = abs ( uex[i] - U[i] )
print ( ' %4d %14.6g %14.6g %14.6g' % ( i, U[i], uex[i], err[i] ) )
plt.plot( x, U, 'bo-',lw=1,label="Numerical ")
plt.plot( xp, up, 'r.',lw=6,label="Real")
plt.savefig ( 'fem1d.png' )
plt.xlabel('$x$', fontsize=14,color='blue')
plt.ylabel('$F(x)$', fontsize=14,color='blue')
plt.legend(fontsize=14)
plt.show ( )
print ( '<--------------------------------------------------->' )
print('使用误差证明该方法是二阶的即 O(h^2)')
print('在无穷范数情况下：')
norm1=max(abs(err))
print('无穷范数下的误差为：%14.6g' % (norm1))
print ( '<--------------------------------------------------->' )
print('在1范数情况下：')
norm2=sum(abs(err))*h
print('1范数下的误差为：%14.6g' % (norm2))
print ( '<--------------------------------------------------->' )
print('在2范数情况下：')
norm3=np.sqrt(sum(abs(err)*(abs(err)).T)*h)
print('2范数下的误差为：%14.6g' % (norm3))
print ( 'PDE——有限差分：中心差商求解：：1D 狄利克雷边值 问题' )
print ( ' Normal end of execution.' )
a,b,n= map(int,input('请输入左端点a,右端点b,单元格个数m：等3个系数：各数字之间空格键隔开\n').split())
start = time.clock()
Fdm_cd_1d(a,b,n)
end = time.clock()
print('计算程序运行的时间\n')
print(end-start)
print('')
Node 数值解 解析解 Error
0 -0 -0 0
1 -0.00315133 -0.00315107 2.54071e-07
2 -0.00629961 -0.0062991 5.07897e-07
3 -0.00944179 -0.00944103 7.61231e-07
4 -0.0125748 -0.0125738 1.01383e-06
5 -0.0156957 -0.0156945 1.26545e-06
6 -0.0188014 -0.0187999 1.51584e-06
7 -0.0218889 -0.0218872 1.76476e-06
8 -0.0249553 -0.0249533 2.01198e-06
9 -0.0279975 -0.0279952 2.25726e-06
10 -0.0310126 -0.0310101 2.50034e-06
11 -0.0339977 -0.0339949 2.74101e-06
12 -0.0369499 -0.0369469 2.97903e-06
13 -0.0398664 -0.0398632 3.21417e-06
14 -0.0427443 -0.0427408 3.44619e-06
15 -0.0455808 -0.0455771 3.67489e-06
16 -0.0483733 -0.0483694 3.90002e-06
17 -0.0511189 -0.0511148 4.12139e-06
18 -0.0538151 -0.0538108 4.33876e-06
19 -0.0564593 -0.0564547 4.55194e-06
20 -0.0590488 -0.059044 4.76072e-06
21 -0.0615812 -0.0615762 4.96489e-06
22 -0.064054 -0.0640488 5.16426e-06
23 -0.0664648 -0.0664595 5.35863e-06
24 -0.0688114 -0.0688058 5.54781e-06
25 -0.0710913 -0.0710856 5.73163e-06
26 -0.0733025 -0.0732966 5.90991e-06
27 -0.0754428 -0.0754367 6.08247e-06
28 -0.0775101 -0.0775039 6.24914e-06
29 -0.0795024 -0.079496 6.40977e-06
30 -0.0814178 -0.0814113 6.56419e-06
31 -0.0832545 -0.0832477 6.71227e-06
32 -0.0850105 -0.0850037 6.85385e-06
33 -0.0866844 -0.0866774 6.9888e-06
34 -0.0882744 -0.0882673 7.11699e-06
35 -0.089779 -0.0897717 7.2383e-06
36 -0.0911967 -0.0911893 7.3526e-06
37 -0.0925262 -0.0925187 7.45979e-06
38 -0.0937662 -0.0937586 7.55976e-06
39 -0.0949155 -0.0949078 7.65242e-06
40 -0.0959729 -0.0959652 7.73767e-06
41 -0.0969375 -0.0969297 7.81544e-06
42 -0.0978083 -0.0978004 7.88565e-06
43 -0.0985845 -0.0985766 7.94823e-06
44 -0.0992654 -0.0992574 8.00312e-06
45 -0.0998501 -0.0998421 8.05027e-06
46 -0.100338 -0.10033 8.08963e-06
47 -0.100729 -0.100721 8.12116e-06
48 -0.101023 -0.101015 8.14484e-06
49 -0.101219 -0.101211 8.16064e-06
50 -0.101317 -0.101309 8.16854e-06
51 -0.101317 -0.101309 8.16854e-06
52 -0.101219 -0.101211 8.16064e-06
53 -0.101023 -0.101015 8.14484e-06
54 -0.100729 -0.100721 8.12116e-06
55 -0.100338 -0.10033 8.08963e-06
56 -0.0998501 -0.0998421 8.05027e-06
57 -0.0992654 -0.0992574 8.00312e-06
58 -0.0985845 -0.0985766 7.94823e-06
59 -0.0978083 -0.0978004 7.88565e-06
60 -0.0969375 -0.0969297 7.81544e-06
61 -0.0959729 -0.0959652 7.73767e-06
62 -0.0949155 -0.0949078 7.65242e-06
63 -0.0937662 -0.0937586 7.55976e-06
64 -0.0925262 -0.0925187 7.45979e-06
65 -0.0911967 -0.0911893 7.3526e-06
66 -0.089779 -0.0897717 7.2383e-06
67 -0.0882744 -0.0882673 7.11699e-06
68 -0.0866844 -0.0866774 6.9888e-06
69 -0.0850105 -0.0850037 6.85385e-06
70 -0.0832545 -0.0832477 6.71227e-06
71 -0.0814178 -0.0814113 6.56419e-06
72 -0.0795024 -0.079496 6.40977e-06
73 -0.0775101 -0.0775039 6.24914e-06
74 -0.0754428 -0.0754367 6.08247e-06
75 -0.0733025 -0.0732966 5.90991e-06
76 -0.0710913 -0.0710856 5.73163e-06
77 -0.0688114 -0.0688058 5.54781e-06
78 -0.0664648 -0.0664595 5.35863e-06
79 -0.064054 -0.0640488 5.16426e-06
80 -0.0615812 -0.0615762 4.96489e-06
81 -0.0590488 -0.059044 4.76072e-06
82 -0.0564593 -0.0564547 4.55194e-06
83 -0.0538151 -0.0538108 4.33876e-06
84 -0.0511189 -0.0511148 4.12139e-06
85 -0.0483733 -0.0483694 3.90002e-06
86 -0.0455808 -0.0455771 3.67489e-06
87 -0.0427443 -0.0427408 3.44619e-06
88 -0.0398664 -0.0398632 3.21417e-06
89 -0.0369499 -0.0369469 2.97903e-06
90 -0.0339977 -0.0339949 2.74101e-06
91 -0.0310126 -0.0310101 2.50034e-06
92 -0.0279975 -0.0279952 2.25726e-06
93 -0.0249553 -0.0249533 2.01198e-06
94 -0.0218889 -0.0218872 1.76476e-06
95 -0.0188014 -0.0187999 1.51584e-06
96 -0.0156957 -0.0156945 1.26545e-06
97 -0.0125748 -0.0125738 1.01383e-06
98 -0.00944179 -0.00944103 7.61231e-07
99 -0.00629961 -0.0062991 5.07897e-07
100 -0.00315133 -0.00315107 2.54071e-07
101 -1.24083e-17 -1.24083e-17 0
<--------------------------------------------------->
使用误差证明该方法是二阶的即 O(h^2)
在无穷范数情况下：
无穷范数下的误差为： 8.16854e-06
<--------------------------------------------------->
在1范数情况下：
1范数下的误差为： 5.20046e-06
<--------------------------------------------------->
在2范数情况下：
2范数下的误差为： 5.77673e-06
PDE——有限差分：中心差商求解：：1D 狄利克雷边值 问题
Normal end of execution.
计算程序运行的时间
0.34405040000001463
3212
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
