使用Julia辛求解器KahanLi8()求解ODE时出现'type Array has no field x'错误

在Julia环境中使用Tsit5()求解器求解常微分方程（ODE）可正常得到结果，但切换为辛求解器KahanLi8()时出现错误：type Array has no field x。完整代码如下：

using LinearAlgebra 
using PyPlot
using Random
using DifferentialEquations

 N = 100
 σ, a, J, K = 10.0, 1.0, 1.0,-0.1;
 ω = zeros(N);  

 function heaviside(x::Real)
 return x >= 0.0 ? 1.0 : 0.0
 end 

 function rhs(du,u,p,t)

  u1 = @view u[1:N]
  du1 = @view du[1:N]
  u2 = @view u[N+1:2*N]
  du2 = @view du[N+1:2*N]
 
  σ,a,J,K=p

  for i in 1:N
    du1[i]=(1/N)*sum(j->( (u1[j]-u1[i])*(1+J*cos(u2[j]-u2[i]))-sign(u1[j]-u1[i]) ),1:N)
    
    du2[i]=ω[i]+ (K/N)* sum(j->( sin(u2[j]-u2[i])*(1-(u1[j]-u1[i])^2/σ^2)*heaviside(σ - 
 abs((u1[j]-u1[i])))),1:N)
end

return du 

end

tspan = (0.0, 200)  # Assuming you want to integrate from 0 to T
dts=0.1 #savetime
p=[σ,a,J,K]
Random.seed!(123)
u0= [(rand() * 20.0 - 10.0, rand() * (2.0 * π) - π) for j in 1:N]
#println(u0)
u0=vcat([x[1] for x in u0], [x[2] for x in u0]);

 prob = ODEProblem(rhs, u0, tspan,p);

 @time sol = solve(prob,KahanLi8(),dt=0.1,saveat=dts);

问题原因 #

辛求解器（如KahanLi8）专门针对可分离哈密顿系统设计，要求状态变量必须是可分离的结构（通常是包含位置q和动量p的数组对，或带有.x字段的复合类型）。当前代码将所有状态变量打包为单一一维数组，不符合辛求解器的输入格式要求，因此触发找不到.x字段的错误。

解决方法 #

需要将状态变量拆分为对应哈密顿系统的q（位置）和p（动量）两部分，调整初始条件和右端函数的格式：

1. 拆分初始条件为两个独立数组，组成元组
2. 修改右端函数，适配拆分后的状态变量输入，分别计算q和p的导数

修改后的完整代码：

using LinearAlgebra 
using PyPlot
using Random
using DifferentialEquations

N = 100
σ, a, J, K = 10.0, 1.0, 1.0,-0.1;
ω = zeros(N);  

function heaviside(x::Real)
    return x >= 0.0 ? 1.0 : 0.0
end 

# 适配辛求解器的右端函数，输入为(q,p)元组
function rhs(du, u, p, t)
    q, p = u  # 拆分状态变量为位置q和动量p
    dq, dp = du  # 拆分导数数组为dq和dp
    
    σ,a,J,K=p

    for i in 1:N
        dq[i]=(1/N)*sum(j->( (q[j]-q[i])*(1+J*cos(p[j]-p[i]))-sign(q[j]-q[i]) ),1:N)
        
        dp[i]=ω[i]+ (K/N)* sum(j->( sin(p[j]-p[i])*(1-(q[j]-q[i])^2/σ^2)*heaviside(σ - abs((q[j]-q[i])))),1:N)
    end

    return du 
end

tspan = (0.0, 200)  
dts=0.1 
p=[σ,a,J,K]
Random.seed!(123)
# 拆分初始条件为q和p两个数组，组成元组
u0_q = [rand() * 20.0 - 10.0 for j in 1:N]
u0_p = [rand() * (2.0 * π) - π for j in 1:N]
u0 = (u0_q, u0_p);

prob = ODEProblem(rhs, u0, tspan,p);

@time sol = solve(prob,KahanLi8(),dt=0.1,saveat=dts);

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