# Julia矩阵计算速度慢的问题

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;
caxis([0,1])
end
``````

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

>> TuringPattern02

>> TuringPattern02

``````

``````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
``````

1 个赞
2 个赞

1 个赞