如何通过表格数据插值得到ε(p)方程？代码报错求助

问题背景 #

需要对能量（ε）和压力（p）的数据集进行插值，得到ε(p)形式的表达式，使用scipy.interpolate.interp1d实现时触发类型错误，无法完成数值计算。

原代码 #

import numpy as np
import math
from scipy.interpolate import interp1d

data_apr = np.loadtxt(r"C:\Users\Ramos\PycharmProjects\pythonProject\\apr.dat")

m_apr = [] 
p_apr = [] #pressure
ε_apr = [] #energy

for i in range(1001):
    m_apr.append(data_apr[i][0])
    p_apr.append(data_apr[i][1])
    ε_apr.append(data_apr[i][2])

εp= interp1d(p_apr, ε_apr)

def f(u, kappa, εp, r):
    λ, v, p = u
    f = np.array([(1-math.exp(u[0])+math.exp(u[0])*(r**2)*kappa*εp)/r, (-1+math.exp(u[0])+math.exp(u[0])*(r**2)*kappa*u[2])/r, -(1/2)*(u[3]+εp)*((-1+math.exp(u[0])+math.exp(u[0])*(r**2)*kappa*u[2])/r)])
    return f

# Runge-Kutta Method

def RK2(u, f, dr, *args):
    u_i = u + dr * f(u, *args) / 2
    u_novo = u + dr * f(u_i, * args)
    return u_novo

#Constants
kappa = 1

dr=0.01
ri = 0.00001
rf = 2
N = int(rf/dr)

u = np.empty((N+1, 3))

#Initial Conditions
λ0 = 1
v0 = 1
p0 = p_apr[0]
u[0] = np.array([λ0, v0, p0])

r = np.arange(0.0001, 1.5, 0.01)

for i in range(N-1):
    u[i+1] = RK2(u[i], f, dr, kappa, εp, r[i])
    if u[i+1][2] < 0 or math.isnan(u[i+1][2]):
        raio = r[i]
        index = i
        break

print("r[index]:", raio)
print("index", index)

报错信息 #

f = np.array([(1-math.exp(u[0])+math.exp(u[0])*(r**2)*kappa*εp)/r, (-1+math.exp(u[0])+math.exp(u[0])*(r**2)*kappa*u[2])/r, -(1/2)*(u[3]+εp)*((-1+math.exp(u[0])+math.exp(u[0])*(r**2)*kappa*u[2])/r)])
TypeError: unsupported operand type(s) for *: 'float' and 'interp1d' 

错误原因与修正方案 #

核心错误 #

interp1d返回的εp是插值函数对象，不是数值，直接将它与浮点数（如r**2、kappa）进行运算会触发类型错误。必须传入当前的压力值p（即u[2]）调用εp(p)，才能得到对应的能量数值。

额外错误 #

原代码中u[3]是无效索引：u是3维数组（存储λ、v、p），索引范围是0-2，此处属于笔误，结合物理意义推测应为能量εp(p)与压力p的和。

修正后的代码 #

import numpy as np
import math
from scipy.interpolate import interp1d

data_apr = np.loadtxt(r"C:\Users\Ramos\PycharmProjects\pythonProject\\apr.dat")

m_apr = [] 
p_apr = [] #pressure
ε_apr = [] #energy

for i in range(1001):
    m_apr.append(data_apr[i][0])
    p_apr.append(data_apr[i][1])
    ε_apr.append(data_apr[i][2])

# 插值函数：输入压力p，返回对应能量ε
εp= interp1d(p_apr, ε_apr)

def f(u, kappa, εp, r):
    λ, v, p = u
    # 调用插值函数获取当前压力对应的能量值
    ε_val = εp(p)
    # 修正u[3]的索引错误，替换为能量+压力
    f = np.array([
        (1 - math.exp(λ) + math.exp(λ) * (r**2) * kappa * ε_val) / r,
        (-1 + math.exp(λ) + math.exp(λ) * (r**2) * kappa * p) / r,
        -(1/2) * (ε_val + p) * ((-1 + math.exp(λ) + math.exp(λ) * (r**2) * kappa * p) / r)
    ])
    return f

# Runge-Kutta Method
def RK2(u, f, dr, *args):
    u_i = u + dr * f(u, *args) / 2
    u_novo = u + dr * f(u_i, * args)
    return u_novo

#Constants
kappa = 1

dr=0.01
ri = 0.00001
rf = 2
N = int(rf/dr)

u = np.empty((N+1, 3))

#Initial Conditions
λ0 = 1
v0 = 1
p0 = p_apr[0]
u[0] = np.array([λ0, v0, p0])

# 确保r数组长度与循环次数匹配
r = np.linspace(ri, rf, N)

for i in range(N-1):
    u[i+1] = RK2(u[i], f, dr, kappa, εp, r[i])
    if u[i+1][2] < 0 or math.isnan(u[i+1][2]):
        raio = r[i]
        index = i
        break

print("r[index]:", raio)
print("index", index)

关键修改点 #

1. 插值函数调用：新增ε_val = εp(p)，用当前压力值p获取对应的能量数值，再参与后续运算。
2. 索引错误修正：将u[3]替换为ε_val + p，符合物理方程中能量密度加压力的常见形式。
3. r数组一致性：将r = np.arange(0.0001, 1.5, 0.01)改为np.linspace(ri, rf, N)，确保r数组长度与Runge-Kutta循环次数匹配，避免索引越界。

内容的提问来源于stack exchange，提问作者Ramos