Python程序中温度更新出现振荡问题的分析和解决方案

在处理温度更新出现振荡问题时，可以考虑以下分析和解决方案：检查温度更新算法是否正确，可能存在错误导致振荡。检查温度更新的步长（时间步长）是否合适，步长过大可能导致振荡。检查系统动力学模型是否准确，可能存在模型不准确导致振荡。

在这里插入图片描述

1、问题背景

在 Python 程序中，通过 class 方法 “update()” 来模拟温度变化时，当 warp 值设置为较高数值（如 1000）时，温度会出现剧烈的振荡。所有温度均以开尔文表示，以简化处理。

2、解决方案

1. 问题分析

问题根源在于积分方案不合理。积分时间步长是 self.warp 个时间单位，即 self.warp 秒，这并不是正确的方法。
需要将方程转换为无量纲形式，即用某种计算单位表示每个项。例如，斯特藩-玻尔兹曼常数和 self.power 可以用这样的单位来衡量，该常数为 1。
然后，应确定对象的特征时间，例如温度达到某种程度的平衡温度所需的时间。如果存在许多这样的对象，则应找到所有特征时间中最小的一个，并将其用作时间的计量单位。
然后，积分时间步长应大约比特征时间小一个数量级，否则将完全错过微分方程的正确解，并最终出现剧烈的振荡。

2. 改进方案

将方程转换为无量纲形式，即用某种计算单位表示每个项。
确定对象的特征时间，例如温度达到某种程度的平衡温度所需的时间。
将积分时间步长设置为大约比特征时间小一个数量级。

3. 代码示例

import math#Critical to the Stefan-Boltzmann equation. Otherwise known as Sigma
BOLTZMANN_CONSTANT = 5.67e-8class GeneratorObject(object):"""Create a new object to run thermal simulation on."""def __init__(self, mass, emissivity, surfaceArea, material, temp=0, power=5000, warp=1):self.tK = temp                                                  #Temperature of the object.     self.mass = mass                                                #Mass of the object.self.emissivity = emissivity                                    #Emissivity of the object. Always between 0 and 1.self.surfaceArea = surfaceArea                                  #Emissive surface area of the object.self.material = material                                        #Store the material name for some reason.self.specificHeat = (0.45*1000)*self.mass                       #Get the specific heat of the object in J/kg (Iron: 0.45*1000=450J/kg)self.power = power                                              #Joules/Second (Watts) input. This is for heating the object.self.warp = warp                                                #Warp Multiplier. This pertains to how KSP's warp multiplier works.def update(self):"""Update the object's temperature according to it's properties."""#This method updates the object's temperature according to heat losses and other factors.# Convert to dimensionless formdimensionless_time = self.warp * self.characteristic_timedimensionless_temperature = self.tK / self.equilibrium_temperature# Update the dimensionless temperaturedimensionless_temperature -= ((self.emissivity * BOLTZMANN_CONSTANT * self.surfaceArea * (math.pow(dimensionless_temperature,4) - math.pow(30+273.15,4))) / self.specificHeat) - (self.power / self.specificHeat)) * dimensionless_time# Convert back to physical unitsself.tK = dimensionless_temperature * self.equilibrium_temperature