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

H.E.G.

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

import math
import decimal as dec
import numpy as np
import scipy as sp

class Model:
    def __dec(self, x):#创建十进制对象
        return dec.Decimal(x)
        
    def __init__(self, a=1.496e11, e=1/60, T=365.25964*sp.constants.day, phi=(66+34/60)/180*math.pi):#初始化
        self.a = self.__dec(a)
        self.e = self.__dec(e)
        self.T = self.__dec(T)
        self.phi = self.__dec(phi)
        self.pi = self.__dec(math.pi)
        
    def Fradium(self, theta):#焦半径计算
        Fradium = self.a * (1 - self.e ** 2) / (1 + self.e * self.__dec(math.cos(theta)))
        return float(Fradium)
        
    def theta_of_t(self, t):#theta(t)的计算，需要用到数值积分，这也是问题的主要出处
        t = self.__dec(t)
        func_of_theta = lambda theta:(1 / (1 + self.e * self.__dec(math.cos(theta)))) ** 2#待积分函数
        res_of_integration = lambda theta:sp.integrate.quad(func_of_theta(theta), 0, theta)#积分结果
        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))#方程
        root = sp.optimize.root(equation, self.pi)#解方程
        try:
            return float(root)
        except:
            print('error')
            return

    def v_of_planet(self, theta):#行星运行速度（如果没算错）
        theta = self.__dec(theta)
        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)#计算式
        return float(v)
        
    def Phi(self, theta):#阳光直射纬度（如果没算错）
        theta = self.__dec(theta)
        Phi = -math.asin(math.cos(theta) * math.cos(self.phi))
        return float(Phi)
        
    def pos_in_xOy(self, theta):#对于的平面直角坐标系中的坐标
        theta = self.__dec(theta)
        Fradium = self.__dec(self.Fradium(theta))
        x = Fradium * math.cos(theta) + self.a * self.e
        y = Fradium * math.sin(theta)
        return (float(x), float(y))

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

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

那么运行结果如下

Warning (from warnings module):
  File "G:\tst\JustACommonModel.py", line 23
    func_of_theta = lambda theta:(1 / (1 + self.e * self.__dec(math.cos(theta)))) ** 2
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.)
Traceback (most recent call last):
  File "G:\tst\JustACommonModel.py", line 54, in <module>
    print(m.theta_of_t(1))
  File "G:\tst\JustACommonModel.py", line 26, in theta_of_t
    root = sp.optimize.root(equation, self.pi)
  ...
ValueError: invalid callable given

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

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

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

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

FishC

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

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

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

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

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

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

def theta_of_t(self, t):
    t = self.__dec(t)
    func_of_theta = lambda theta: (self.__dec(1) / (self.__dec(1) + self.e * self.__dec(math.cos(theta)))) ** 2
    res_of_integration = sp.integrate.quad(func_of_theta, 0, theta)
    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))
    root = sp.optimize.root(equation, self.pi)
    try:
        return float(root.x[0])
    except:
        print('error')
        return None

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

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

歌者文明清理员

nb

liuhongrun2022

给我整不会了