Julia耦合ODE并行求解优化咨询（M1 Mac环境）

问题背景 #

我正在求解耦合ODE，仅展示问题的一小部分。由于网格规模过大，常规方法无法运行，因此尝试在Julia中实现并行化，但在搭载M1芯片、8GB内存的Mac环境中运行代码时始终崩溃。需求是求解整个u和v网格上的R值：先通过牛顿法求解u=-100、v=0处的R_0，再用RK4格式求解其余网格点，烦请帮忙排查代码问题。

代码片段 #

function f(R, v)
    numerator = R * exp(1 / (3 - 50 * R))
    denominator = (2 + (200 * R) / 3)^(4 / 9) * (2 - (100 * R) / 3)^(5 / 9)
    return v - (numerator / denominator)
end


function df(R, v)
    # Approximate the derivative using central difference method
    h = 1e-6
    return (f(R + h, v) - f(R - h, v)) / (2 * h)
end

# Newton's method implementation
function newton(f, df, R0, v; tol=1e-6, maxiter=1000)
    R = R0
    for i in 1:maxiter
        R_new = R - f(R, v) / df(R, v)
        if abs(R_new - R) < tol
            return R_new
        end
        R = R_new
    end
    error("Newton's method did not converge")
end

# Array to store final converged values of R
# This is the value along the intial hypersurface
R_0 = Float64[] 
v = LinRange(0.0000, 150, 110001)
# Initial conditions
let R_initial = 0.025  # Initial guess for R
    for v in v
    R_solution = newton(f, df, R_initial, v)
    R_initial = R_solution
    push!(R_0, R_solution) # Update R0 to the last converged value of R
    end
end
using PlotlyJS
plot(R_0, v)

function Initial_Function(R)
    H_0 = R^2 * (0.0525 - R)^2
end

H_initial = Float64[]
for R in R_0
    if R <= 0.0525
        H = Initial_Function(R)
        push!(H_initial, H)
    end
    if R > 0.0525
        H = 0
        push!(H_initial, H)
    end
end
using PlotlyJS
plot(R_0, H_initial)
# H, R, v are defined
# We run the RK4 scheme to calculate the value of R across the grid.

u = range(-100, 10, length = 110001)
v_range = range(0,150,length = 110001)
using Distributed
addprocs(4)
@everywhere using SharedArrays
R_grid = SharedArray{Float64}(110001,110001)
R_grid[1,:] = R_0

function dR_du(u, R)
    g = (R^2/6)*(u*R + (u^2 * R^2)/9)
    return (- 0.5 * g)
end

function rk4(u, R, du)
    k1 = dR_du(u, R)
    k2 = dR_du(u + 0.5*du, R+0.5*du*k1)
    k3 = dR_du(u + 0.5*du, R + 0.5*du*k2)
    k4 = dR_du(u + du, R+ du*k3)
    return R + (du/6) * (k1 + 2*k2 + 2*k3 + k4)
end
Threads.@threads for j in 1:110001 
    for i in 2:110001
        u_vals = u[i-1]
        du = u[i] - u_vals
        R_grid[j, i] = rk4(u_vals, R_grid[j, i-1], du)
    end
end

问题排查点 #

• 内存占用超限：定义的R_grid是110001×110001的Float64数组，单个Float64占8字节，总内存需求约为96GB，远超过8GB的机器内存，这是崩溃的核心原因。即使使用SharedArray，总内存需求不变，只是分配到不同进程，你的机器根本无法承载。
• 并行模式混用：同时使用Distributed的addprocs和Threads.@threads，两种并行模式（分布式多进程vs多线程）混用可能引发资源冲突，进一步加剧内存问题。
• 数值稳定性隐患：Newton法中用中心差分近似导数，当R接近0.06时（分母中2 - 100R/3趋近于0），f(R,v)会出现数值奇异，可能导致Newton法发散或导数计算溢出，后续RK4迭代也会受影响。
• 冗余导入：代码中两次调用using PlotlyJS，属于冗余操作，虽不直接导致崩溃，但影响代码整洁。

优化建议 #

• 缩小网格规模：先使用小网格（如1001×1001）验证逻辑正确性，再根据内存情况逐步调整；或采用分块计算、稀疏存储，只处理当前需要的部分数据，避免一次性加载整个大网格。
• 选择单一并行模式：优先用多线程（Threads.@threads）处理CPU密集型任务，无需addprocs和SharedArray；若必须分布式，确保每个进程只处理部分数据块，降低单进程内存占用。
• 替换数值导数为解析导数：手动推导f(R,v)的解析导数，替换中心差分，既提升计算效率，又避免奇异点的导数计算误差。
• 预分配数组：初始化数组时直接预分配空间（如R_0 = Vector{Float64}(undef, length(v))），比push!更高效，减少内存碎片。

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