【发布时间】：2012-10-09 01:26:32
【问题描述】：
我正在寻求有关我在 Mathematica 中编写的 Poisson Solver 的帮助。插入数组的代码很长，但完整的细节可以在http://pastebin.com/uSrSDcW6找到。
我正在使用源自泊松方程的中心差分法计算给定电荷密度的电压。计算电压后，我测试数据集的收敛性。我将收敛阈值设置为 10^-1000+ 的数量级。我将循环设置为在 10000 次迭代后退出，以防万一出现问题，作为故障保险。我有一个循环计数器以保持理智。只要收敛阈值设置为 10^-100，程序似乎运行良好。
我的问题是：无论我更新阈值，例如 10^-100、10^-150，计算在 633 次迭代后停止并退出循环。我将不胜感激，我完全被困住了。我已将 cmets 添加到程序中，该程序应该对本论坛上的任何人都具有解释性。再说一次，我知道这个描述是有限的，所以请查看附件的 url http://pastebin.com/uSrSDcW6 以获得完整的程序。
*更新10/9/12***我已将问题隔离到 16 位机器精度。我需要将它打开到我的机器最大精度为 10^309。 Mathematica 帮助很少说明如何做到这一点。例如“N [机器精度，50]”。我将在我的程序中将其设置在哪里以将其应用于所有计算？如果有帮助，我会在此处粘贴循环
Vnew / Vold / RHO 是 10x10x34 矩阵
Epsilon 是一个常数
（将 ConvergenceLoop 初始化为 O - 如果需要，这将作为故障安全退出循环）
收敛环 = 0;
(将收敛初始化为零)
收敛 = 0;
而[Convergence == 0 && ConvergenceLoop
(遍历所有 i,j,k 元素，计算新的电压值)
Do[Vnew[[i]][[j]][[k]] = (1/(2/deltaX^2 + 2/deltaY^2 + 2/deltaZ^2)) *(((体积[[i + 1]][[j]][[k]] +
Vold[[i - 1]][[j]][[k]])/(deltaX^2)) + ((Vold[[i]][[j + 1]][[k]] +
Vold[[i]][[j - 1]][[k]])/(deltaY^2)) + ((Vold[[i]][[j]][[k + 1]] +
Vold[[i]][[j]][[k - 1]])/(deltaZ^2)) + ((Rho[[i]][[j]][[k]]/Epsilon))), {i, 2, 9}, {j, 2,9}, {k, 2, 33}];
(假设收敛，因此当测试达到超过定义的收敛阈值的第一个值时触发循环)
收敛 = 1;
(这是收敛测试。用户定义的收敛阈值)
Do[If[Vold[[i]][[j]][[k]] == 0, Null,
If[(Vnew[[i]][[j]][[k]] - Vold[[i]][[j]][[k]])/Vold[[i]][[j]][[k]] > .0000001, Convergence = 0;
(*This is purely diagnostic. I added a Tracker point to follow the convergence along.
用户在任何元素上定义*)
If[i == 5 && j == 5 && k == 10,
Print[ "Tracker Point" (Vnew[[i]][[j]][[k]] -
Vold[[i]][[j]][[k]])/Vold[[i]][[j]][[k]]], Null],Null]], {i, 2, 9}, {j, 2, 9}, {k, 2, 33}];
(忽略第一次迭代，直到 Vnew 和 Vold 不为零)
如果[ConvergenceLoop
（迫使 Vold 与 Vnew 一起进化）
体积 = Vnew;
收敛循环 ++;]
（为未来的规划目的添加了 SessionTime）
如果[ConvergenceLoop == 10000，
Print["达到收敛循环限制。" (SessionTime[]/3600)],
Print["未达到收敛循环限制。"]];
(我们打破了while循环，意味着我们的数据已经收敛，所以打印收敛的值)
如果[收敛 == 1，
打印[ConvergenceLoop“恭喜你收敛了！” MatrixForm [Vnew]], Print["没有收敛！"]];
【问题讨论】：
-
我没心情去挖掘你的代码，但也许你需要使用任意精度，而你正在使用机器精度？如果您正在比较两个应该不同但不是由于舍入误差而导致的数字，则可能会出现这种情况。
-
了解 - 我已将问题隔离到 16 位机器精度。我需要将它打开到我的机器最大精度为 10^309。 Mathematica 帮助很少说明如何做到这一点。例如“N [机器精度，50]”。我应该在我的程序中的哪里设置它以将其应用于所有计算？
标签： wolfram-mathematica solver poisson