各位好！求助一个PYTHON做偏微分数值计算精度的问题

wuxianke

在网上找了个算例，解二维抛物线偏微分方程，按算法跑了1600步，结果精确值为0的地方变成0.002几。但我看算例上用MATLAB做，计算值才10^-16。为什么差别这么大呢，是我的代码有问题还是，需要用什么方法提高PYTHON的精度呢？

希望有大神可以解惑！

附上自己的代码和算例链接

1. # -*- coding: utf-8 -*-
   
2. """
   
3. Created on Fri Nov  8 18:04:41 2019
   
4. 
   
5. @author: Administrator
   
6. """
   
7. from __future__ import division
   
8. import numpy as np
   
9. import math
   
10. import matplotlib.pyplot as plt
    
11. 
    
12. 
    
13. 
    
14. #计算初始U(j,k)
    
15. 
    
16. h = 1/40  #空间步长
    
17. t = 1/1600 #空间步长
    
18.     #a2 = 1/16 #系数
    
19. miu = t/h**2
    
20. F1 = np.zeros((38,38))
    
21. F2 = np.zeros((40,40))
    
22. def doit():
    
23. 
    
24.     a1,b1,c1 = construct_step1_A()
    
25.     a2,b2,c2 = construct_step2_A()
    
26.     U0 = initialize()
    
27.     T = 0
    
28.     er = []
    
29.     while T<30:
    
30.         Uhalf = forward_half_1(a1,b1,c1,U0)
    
31.         U1 = forward_half_2(a2,b2,c2,Uhalf)
    
32.         U0 = U1
    
33.         T += 1
    
34. 
    
35.     return U0
    
36. 
    
37. #计算精确值
    
38. def realsolve(j,k,t):
    
39.     return math.sin(math.pi*j*h)*math.cos(math.pi*k*h)*math.e**(-math.pi**2/8*t)
    
40. 
    
41. def initialize():
    
42.     U = np.zeros((40,40))
    
43.     for i in range(40):
    
44.         for j in range(40):
    
45.             U[i,j] = realsolve(i,j,0)
    
46.         
    
47.     return U
    
48. 
    
49. 
    
50. #追赶法求解
    
51. def thomas_solve(a,b,c,f):
    
52.     x = np.zeros_like(f)
    
53.     y = np.zeros_like(f)   
    
54.     u = np.zeros_like(f)   
    
55.     l = np.zeros((len(f)-1,1))   
    
56.     u[0] = b[0]
    
57.     for i in range(1,len(f)):
    
58.         l[i-1] = a[i-1]/u[i-1]
    
59.         u[i] = b[i] - l[i-1]*c[i-1]
    
60.    
    
61.     y[0] = f[0]
    
62.     for i in range(1,len(f)):
    
63.         y[i] = f[i] - l[i-1] * y[i-1]
    
64.    
    
65.     x[len(x)-1] = y[len(x)-1]/u[len(x)-1]
    
66.     for i in range(len(x)-2,-1,-1):
    
67.         x[i] = (y[i] - c[i]*x[i+1])/u[i]
    
68.    
    
69.     return x
    
70.       
    
71.    
    
72. #初始化第一步A矩阵
    
73. def construct_step1_A():   
    
74.     a = np.ones((37,1))*(-1/32*miu)
    
75.     c = np.ones((37,1))*(-1/32*miu)
    
76.     b = np.ones((38,1))*(1+1/16*miu)
    
77.    
    
78.     return a,b,c
    
79.    
    
80. def construct_step2_A():
    
81.     a = np.ones((39,1))*(-1/32*miu)
    
82.     c = np.ones((39,1))*(-1/32*miu)
    
83.     b = np.ones((40,1))*(1+1/16*miu)
    
84.    
    
85.     b[0] = 1
    
86.     b[39] = 1
    
87.     c[0] = -1
    
88.     a[38] = -1
    
89.    
    
90.     return a,b,c
    
91.    
    
92. 
    
93. #初始化0时刻向量矩阵F
    
94. def construct_step1_F(U0):   
    
95.     for j in range(38):
    
96.         for k in range(38):
    
97.             F1[j,k] = miu/32*(U0[j+1,k+1-1] - 2*U0[j+1,k+1] + U0[j+1,k+1+1]) + U0[j+1,k+1]
    
98.    
    
99.     return F1
    
100. 
     
101. def construct_step2_F(U1):   
     
102.     for j in range(1,39):
     
103.         for k in range(1,39):
     
104.             F2[j,k] = miu/32*(U1[j-1,k] - 2*U1[j,k] + U1[j+1,k]) +U1[j,k]
     
105.    
     
106.     return F2
     
107. 
     
108. 
     
109. def forward_half_1(a1,b1,c1,U0):
     
110.     F1 = construct_step1_F(U0)
     
111.     Uhalf = np.zeros_like(U0)
     
112.     for k in range(1,Uhalf.shape[0]-1):
     
113.         Uhalf[1:Uhalf.shape[0]-1,k] = thomas_solve(a1,b1,c1,F1[:,k-1])
     
114. 
     
115.     return Uhalf
     
116. 
     
117. 
     
118. def forward_half_2(a2,b2,c2,Uhalf):
     
119.     F2 = construct_step2_F(Uhalf)
     
120.     U1 = np.zeros_like(Uhalf)
     
121.     for j in range(U1.shape[0]):
     
122.         U1[j,:] = thomas_solve(a2,b2,c2,F2[j,:])
     
123.       
     
124.     return U1
     
125. 
     
126. 
     
127. 
     
128.             
     
129.    
     
130. 
     

复制代码

算例地址

wp231957

这方面MATLAB才是专业吧

wuxianke

wp231957 发表于 2019-11-12 17:51
这方面MATLAB才是专业吧

想用PYTHON做做。。如果不行就去用MATLAB做。。

wuxianke

求帮助！