关于用python建立的某一数学模型的求助

H.E.G.

咳咳，让我算算，3年没泡论坛了，可能出现帖子格式不对，给大家带来麻烦，见谅
回归正题
代码：

1. import math
   
2. import decimal as dec
   
3. import numpy as np
   
4. import scipy as sp
   
5. 
   
6. class Model:
   
7.     def __dec(self, x):#创建十进制对象
   
8.         return dec.Decimal(x)
   
9.         
   
10.     def __init__(self, a=1.496e11, e=1/60, T=365.25964*sp.constants.day, phi=(66+34/60)/180*math.pi):#初始化
    
11.         self.a = self.__dec(a)
    
12.         self.e = self.__dec(e)
    
13.         self.T = self.__dec(T)
    
14.         self.phi = self.__dec(phi)
    
15.         self.pi = self.__dec(math.pi)
    
16.         
    
17.     def Fradium(self, theta):#焦半径计算
    
18.         Fradium = self.a * (1 - self.e ** 2) / (1 + self.e * self.__dec(math.cos(theta)))
    
19.         return float(Fradium)
    
20.         
    
21.     def theta_of_t(self, t):#theta(t)的计算，需要用到数值积分，这也是问题的主要出处
    
22.         t = self.__dec(t)
    
23.         func_of_theta = lambda theta:(1 / (1 + self.e * self.__dec(math.cos(theta)))) ** 2#待积分函数
    
24.         res_of_integration = lambda theta:sp.integrate.quad(func_of_theta(theta), 0, theta)#积分结果
    
25.         equation = lambda theta:float(res_of_integration(theta) - 2 * self.pi * math.sqrt(1 - self.e**2) * t / (self.T * (1 - self.e ** 2) ** 2))#方程
    
26.         root = sp.optimize.root(equation, self.pi)#解方程
    
27.         try:
    
28.             return float(root)
    
29.         except:
    
30.             print('error')
    
31.             return
    
32. 
    
33.     def v_of_planet(self, theta):#行星运行速度（如果没算错）
    
34.         theta = self.__dec(theta)
    
35.         v = 2 * self.pi * self.a * self.__dec(math.sqrt(1 + self.e ** 2 + 2 * self.e * self.__dec(math.cos(theta)))) / (math.sqrt(1 - self.e ** 2) * self.T)#计算式
    
36.         return float(v)
    
37.         
    
38.     def Phi(self, theta):#阳光直射纬度（如果没算错）
    
39.         theta = self.__dec(theta)
    
40.         Phi = -math.asin(math.cos(theta) * math.cos(self.phi))
    
41.         return float(Phi)
    
42.         
    
43.     def pos_in_xOy(self, theta):#对于的平面直角坐标系中的坐标
    
44.         theta = self.__dec(theta)
    
45.         Fradium = self.__dec(self.Fradium(theta))
    
46.         x = Fradium * math.cos(theta) + self.a * self.e
    
47.         y = Fradium * math.sin(theta)
    
48.         return (float(x), float(y))
    

复制代码

对于上面的数学模型有以下几个问题：
1.如果调用下面语句：

1. if __name__=='__main__':
   
2.     m = Model()
   
3.     print(m.theta_of_t(1))

复制代码

那么运行结果如下

1. Warning (from warnings module):
   
2.   File "G:\tst\JustACommonModel.py", line 23
   
3.     func_of_theta = lambda theta:(1 / (1 + self.e * self.__dec(math.cos(theta)))) ** 2
   
4. DeprecationWarning: Conversion of an array with ndim > 0 to a scalar is deprecated, and will error in future. Ensure you extract a single element from your array before performing this operation. (Deprecated NumPy 1.25.)
   
5. Traceback (most recent call last):
   
6.   File "G:\tst\JustACommonModel.py", line 54, in <module>
   
7.     print(m.theta_of_t(1))
   
8.   File "G:\tst\JustACommonModel.py", line 26, in theta_of_t
   
9.     root = sp.optimize.root(equation, self.pi)
   
10.   ...
    
11. ValueError: invalid callable given

复制代码

应该是积分函数的调用之类的问题吧，代码较硬，难以解释，有志解决问题的壮士敬请谅解

2:.代码实在过硬，难以解释，就是上面说的

3.精度不足，而利用decimal中的Decimal()创建实例化对象后，会有类型不统一而精度下降或抛出错误
等等诸多问题

还有待有志之士前来挑战这几个难得让我晕头转向的bug
如果代码还存在除上面提到以外的问题还请指正

FishC

嘿，代码看起来有点复杂，但我会尽力解释并帮助你解决问题。

1. 首先，关于DeprecationWarning的警告，它是由于在调用函数quad时，传入了一个lambda函数而不是一个可调用对象。你可以将这个lambda函数改为func_of_theta本身，这样警告就不会出现了。修改后的代码如下：

1. res_of_integration = sp.integrate.quad(func_of_theta, 0, theta)

复制代码

2. 第二个问题与第一个问题类似，也是关于函数的调用。在使用sp.optimize.root函数时，你应该传入一个可调用对象，而不是一个函数的返回值。因此，将equation函数作为参数传入root函数即可：

1. root = sp.optimize.root(equation, self.pi)

复制代码

3. 关于精度问题，你可以使用decimal库提供的Decimal类型来保持计算精度。在初始化中，你已经将a、e、T和phi转换成了Decimal类型。但是，在theta_of_t函数中，你需要确保所有的中间计算都使用Decimal类型。所以，你可以将func_of_theta和equation函数中的计算都改为使用Decimal类型，以保持精度。修改后的代码如下：

1. def theta_of_t(self, t):
   
2.     t = self.__dec(t)
   
3.     func_of_theta = lambda theta: (self.__dec(1) / (self.__dec(1) + self.e * self.__dec(math.cos(theta)))) ** 2
   
4.     res_of_integration = sp.integrate.quad(func_of_theta, 0, theta)
   
5.     equation = lambda theta: float(res_of_integration[0] - 2 * self.pi * math.sqrt(1 - self.e**2) * t / (self.T * (1 - self.e ** 2) ** 2))
   
6.     root = sp.optimize.root(equation, self.pi)
   
7.     try:
   
8.         return float(root.x[0])
   
9.     except:
   
10.         print('error')
    
11.         return None

复制代码

希望这些修改能够解决你的问题！如果还有其他疑问，请随时提问。

以上回复来自 FishC AI 小助理，如未能正确解答您的问题，请继续追问。

歌者文明清理员

nb

liuhongrun2022

给我整不会了