公告:1)网站程序升级:Q2A升级到1.8.6,Wordpress升级到5.7.2
2)修复了头像加载慢与提交问题反应慢等问题
2021-06-16

欢迎来到 Mathematica 问答社区

提问时请贴上文本代码

语法高亮:在编辑器中点击

被禁止的话题:广告破解

请阅读:《提问的智慧》

备用域名:mma.ooo

支持LaTex数学公式
行内公式标识符:\$ 或“$\backslash ($”+“$\backslash )$”,
行间公式标识符:\$\$ 或 “$\backslash [$”+“$\backslash ]$”

社区建议QQ群:365716997

分类

+1 投票
1.2k 浏览

我的问题是这样的,k[x_, t_] := t - x;
f[x_] := x;
t0 = 0.;
tn = 4.;
n = 400;
h = (tn - t0)/n;
b = f[#] & /@ (h*Range[0, n]);
m = SparseArray[{{i_, j_} /; j < i && j != 1 -> -h*
      k[h*i, h*j], {i_, 1} /; i != 1 -> -0.5*h*k[h*i, 0], {i_, i_} -> 
     1 - 0.5*h*k[h*i, h*i]}, {n + 1, n + 1}];
result1 = LinearSolve[m, b];

分类:矩阵 | 用户: 随堂测验 (606 分)
修改于 用户:随堂测验

1个回答

+4 投票
 
已采纳

有意思的问题。按理来说这个规模的矩阵生成,就算用for循环遍历往里填数也不应该这么慢,也不知道mma是抽了什么风。

其中一个加速的方法是向量化,思路很简单我就直接上代码了

m = DiagonalMatrix[1 - 0.5 h*k[h*#, h*#] &[Range[n + 1]]];
m[[2 ;; -1, 1]] = -0.5*h*k[h*Range[2, n + 1], 0];
Do[m[[j + 1 ;; -1, j]] = -h*k[h*Range[j + 1, n + 1], h*j], {j, 2, 
  n + 1}];
m = SparseArray@m;

效果如下

然而这种方法其实并不能解决问题,因为当n=20000的时候这样生成系数矩阵会爆内存(我16G内存的电脑直接卡死机了)

事实上作为一个下三角矩阵,想要求解的话根本没有必要吧完整的矩阵生成出来,直接一行一行的迭代就好了,代码如下:

ans = Module[{x, alist},
    x = ConstantArray[0., n + 1];
    x[[1]] = b[[1]]/(1 - 0.5 h*k[h, h]);
    Do[
     alist = -h*k[h*i, h*Range[i - 1]];
     alist[[1]] = -0.5*h*k[h*i, h*i];
     x[[i]] = (b[[i]] - alist.x[[1 ;; i - 1]])/(
      1 - 0.5 h*k[i h, i h]), {i, 2, n + 1}];
    x]; // AbsoluteTiming

n=4000时效果如下图,可以看到速度快乐很多,而且内存占用也从O(n^2)降到了O(n)。按照这种做法n=20000时用时也只有3秒多

用户: 无影东瓜 (1.2k 分)
采纳于 用户:随堂测验
非常感谢瓜哥,每次看你的解答都很享受。
“按理来说这个规模的矩阵生成,就算用for循环遍历往里填数也不应该这么慢,也不知道mma是抽了什么风。”并没有抽风,在这个问题上用SparseArray效率不佳我并不意外,因为这个阵并不稀疏。(有约一半的非零元素。)然后,这个问题迭代法无疑是正确的思路。(如果不从根本上去修改求解方法的话。)
确实啊,不懂算法编不好程序啊,瓜哥慧眼直接将空间复杂度从O(n^2)降到了O(n)。
赞,学习啦~嗯,可能速度慢的原因确实是那个矩阵不是稀疏矩阵
...