Julia矩阵计算速度慢的问题

综合讨论区 Julia入门
1

这里分别用Julia和Matlab求解图灵斑问题，两边的代码都是类似的。
Julia的：

using Plots
using LinearAlgebra, SparseArrays

const L = 1.0
const dt = 0.1
const tmax = 100.0
const N = 128
const Dᵤ = 0.4*(L/N)^2         # 为使求解稳定
const Dᵥ = 0.09*(L/N)^2
const F = 0.09
const K = 0.06

function main()
    # 构造一维二阶导算子D₂
    h = L/N
    D₂ = sparse(1/h^2 * diagm(1 => ones(N-1), 0 => -2*ones(N), -1 => ones(N-1)))
    D₂[1, N] = 1/h^2; D₂[N, 1] = 1/h^2     # 周期性边界处理

    # 初值
    u = rand(N, N)
    v = rand(N, N)
    
    # 求解
    kt = length(1:dt:tmax)
    @time for i in 1:kt 
        u = u + dt*(Dᵤ*(u*D₂' + D₂*u) - u.*v.^2 + F*(1 .- u))
        v = v + dt*(Dᵥ*(v*D₂' + D₂*v) + u.*v.^2 - (F + K)*v )   # 直接用解得的新u
    end

    # uv分布绘图
    X = LinRange(0, L, N)
    Y = LinRange(0, L, N)
    fig1 = heatmap(X, Y, u, title="u", color=:lisbon)
    fig2 = heatmap(X, Y, v, title="v")
    fig3 = plot(fig1, fig2,
    layout=(2, 1),
    aspect_ratio=1,
    xlims=(0, L),
    ylims=(0, L),
    clims=(0, 1),
    size=(500,800))
end
main()

Matlab的：

clc
clear
close all
% Gray Scott模型
%% 主函数
main()
function main()
%% 离散设置
L = 1.0;
dt = 0.1;
tmax = 100.0;
N = 128;
%% 常数
Du = 0.4*(L/N)^2;         % 为使求解稳定
Dv = 0.09*(L/N)^2;
F = 0.09;
K = 0.06;
%% 构造一维二阶导算子
h = L/N;
D2d = sparse(1/h^2 * (diag(ones(N-1,1),1) + diag(-2*ones(N,1),0) + diag(ones(N-1,1),-1)));
D2d(1, N) = 1/h^2; D2d(N, 1) = 1/h^2;     % 周期性边界处理
%% 初值
u = rand(N, N);
v = rand(N, N);
%% 求解
kt = length(1:dt:tmax);
tic
for i = 1:kt
    u = u + dt*(Du*(u*D2d' + D2d*u) - u.*v.^2 + F*(1 - u));
    v = v + dt*(Dv*(v*D2d' + D2d*v) + u.*v.^2 - (F + K)*v);
end
toc
%% 绘制u分布图
x = h*(1:N);
[X,Y]=meshgrid(x,x);
pcolor(X,Y,u);
colorbar;
shading interp
caxis([0,1])
end

主要问题出现在注释的"求解"部分，可以看到两边的矩阵计算都是一样的，拉普拉斯算子都是稀疏矩阵格式，但是Julia的求解用时却是Matlab的3倍。
Julia：

  1.212730 seconds (41.62 k allocations: 2.541 GiB, 10.33% gc time)

  1.270115 seconds (41.62 k allocations: 2.541 GiB, 11.55% gc time)

  1.329900 seconds (41.62 k allocations: 2.541 GiB, 11.39% gc time)

Matlab：

>> TuringPattern02
时间已过 0.370999 秒。
>> TuringPattern02
时间已过 0.365901 秒。
>> TuringPattern02
时间已过 0.377799 秒。

不太清楚是什么导致了这样的差距。
本人还是刚刚入门Julia这门语言，还望大佬能答疑解惑，也希望能帮助到有类似困惑的人，谢谢。

2

后来注意到Julia的矩阵转置不是 ’ (共轭转置)，而是transpose()函数，但改了对性能的帮助不明显。
现在我将有限差分求解这段拆解为for循环，不再用矩阵计算了，性能得到很大的提升：

using Plots

const L = 1.0
const dt = 0.1
const tmax = 100.0
const N = 128
const Dᵤ = 0.4*(L/N)^2         # 为使求解稳定
const Dᵥ = 0.09*(L/N)^2
const F = 0.09
const K = 0.06

function main()
    # 初值
    u = rand(N, N)
    v = rand(N, N)
    
    # 求解
    kt = length(1:dt:tmax)
    h = L/N
    @time for k in 1:kt
        for i in 1:N
            for j in 1:N
                u[i, j] = u[i, j] + dt*(Dᵤ*1/h^2*(u[looplc(i-1), j] + u[looplc(i+1), j] + u[i, looplc(j-1)] + u[i, looplc(j+1)] - 4*u[i, j]) -
                        u[i, j]*v[i, j]^2 + F*(1 - u[i, j]))
                v[i, j] = v[i, j] + dt*(Dᵥ*1/h^2*(v[looplc(i-1), j] + v[looplc(i+1), j] + v[i, looplc(j-1)] + v[i, looplc(j+1)] - 4*v[i, j]) +
                        u[i, j]*v[i, j]^2 - (F + K)*v[i, j])
            end
        end
    end

    # uv分布绘图
    X = LinRange(0, L, N)
    Y = LinRange(0, L, N)
    fig1 = heatmap(X, Y, u, title="u", color=:lisbon)
    fig2 = heatmap(X, Y, v, title="v")
    fig3 = plot(fig1, fig2,
    layout=(2, 1),
    aspect_ratio=1,
    xlims=(0, L),
    ylims=(0, L),
    clims=(0, 1),
    size=(500,800))
end

# 循环索引(周期边界)
function looplc(m)
        if m == 0
                N
        elseif m == N+1
                1
        else
                m
        end 
end
main()

3次求解用时：

  0.162914 seconds

  0.258794 seconds

  0.250695 seconds

但是我还是想知道为什么写成矩阵形式就比较慢，而且很多时候也希望将算子化为矩阵形式好向量化编程，希望能有大佬指点一下 :slightly_smiling_face:

1 个赞
3

image
image1384×2092 478 KB

2 个赞
4

以后论坛应该开一个 ChatGPT 自动回复的机器人 :joy:

1 个赞
5

原来是这样，谢谢大佬。不过询问GPT我还真是没想到 :smile:

6

感觉关键可能是 SparseArray？你重写之后用的是 DenseArray 了。