基于松弛法（有限差分）的非线性耦合BVP求解问题咨询

首次在此发帖，因LaTeX无法正常使用，故将部分内容以图片形式发布：

松弛法参考资料为Veseley所著《Computational Physics》。我采用该方法求解给定的非线性耦合ODE，代码核心参数为网格点数n、迭代次数m及松弛参数w。仅采用双侧Dirichlet边界条件时，可轻松达到O(10^{-15})量级的残差容限，但目前遇到以下问题：

• 仅少数n、m、w参数组合可正常运行，调整初始常数（如u(a)）或参数以提升精度时，需重新寻找合适的w，代码稳定性差，易出现病态矩阵；
• 代码正常运行时，Neumann边界条件导致残差容限无法达到Dirichlet条件下的水平，难以同时满足控制方程与边界条件。

以下为可正常运行的Mathematica代码及求解结果图：

(*Equation setup*)
Needs["VariationalMethods`"]
u0 = 10;
m = (1/u[x] - 1/u0);
f = 1 - m z[x]^(d + 1);
L = (Sqrt[-f u'[x]^2 - 2 u'[x] z'[x] + 1]/z[x]^d);(*Lagrangian*)
eq1 = EulerEquations[L, u[x], x];(*Euler-Lagrange equations*)
eq2 = EulerEquations[L, z[x], x];
s = Solve[{eq1, eq2} /. d -> 3, {u''[x], z''[x]}][[1]] // Simplify;
eq01 = (u''[x] - s[[1, 2]]);(*eq01 & eq02 are the nonlinear coupled ODEs*)
eq02 = (z''[x] - s[[2, 2]]);

(*Finite difference and construction of the residual*)
n = 800;(*number of grid points*)
h = (b - a)/(n - 1);(*Step size*)
a = 10^-1;(*a & b are the domain*)
b = 5;
ui = 5;(*ui, zi, are the Dirichlet boundary conditions*)
zi = 10^-1;
uf = 9;(*uf, zf, are the test values for the final grid points u[n], z[n]*)
zf = 5;
up = 0;(*up & zp are the Neumann boundary conditions which is incorporated into the residual down below*)
zp = 0;
rule = {u''[x] -> ((u[i + 1] - 2 u[i] + u[i - 1])/h^2), u'[x] -> ((u[i + 1] - u[i - 1])/(2 h)), u[x] -> u[i], z''[x] -> ((z[i + 1] - 2 z[i] + z[i - 1])/h^2), z'[x] -> ((z[i + 1] - z[i - 1])/(2 h)), z[x] -> z[i]};(*finite difference rule for the interior points*)
resid1 = h^2 eq01 /. rule;(*resid1 & resid2 are the interior points' residuals*)
resid2 = h^2 eq02 /. rule;
residbound1 = (-u[n - 2] + 8 u[n - 1] - 7 u[n] + 6 up h) - 2 h^2 s[[1, 2]] /. {d -> 3, u[x] -> u[n], z[x] -> z[n], u'[x] -> up, z'[x] -> zp};(*residbound1 & residbound2 are the residuals at grid point n*)
residbound2 = (-z[n - 2] + 8 z[n - 1] - 7 z[n] + 6 zp h) - 2 h^2 s[[2, 2]] /. {d -> 3, u[x] -> u[n], z[x] -> z[n], u'[x] -> up, z'[x] -> zp};

(*Construction of sparse matrix*)
parresid1 = {D[resid1, u[i - 1]], D[resid1, z[i - 1]], D[resid1, u[i]], D[resid1, z[i]], D[resid1, u[i + 1]], D[resid1, z[i + 1]]};(*partial derivative of the residual*)
parresid2 = {D[resid2, u[i - 1]], D[resid2, z[i - 1]], D[resid2, u[i]], D[resid2, z[i]], D[resid2, u[i + 1]], D[resid2, z[i + 1]]};
parresidbound1 = D[residbound1, {{u[n - 2], z[n - 2], u[n - 1], z[n - 1], u[n], z[n]}}];(*partial derivative of the boundary residual*)
parresidbound2 = D[residbound2, {{u[n - 2], z[n - 2], u[n - 1], z[n - 1], u[n], z[n]}}];
mat = {parresid1, parresid2};
sparseresidual = Normal[SparseArray[Table[Band[{2 (i - 2) + 1, 2 (i - 2) + 1}] -> {mat}, {i, 2, n - 1}]]];
sparse = Join[{Join[{1}, ConstantArray[0, 2 n - 1]]}, {Join[{0, 1}, ConstantArray[0, 2 n - 2]]}, sparseresidual, {Join[ConstantArray[0, 2 n - 6], parresidbound1]}, {Join[ConstantArray[0, 2 n - 6], parresidbound2]}];(*sparse matrix*)

(*Solve the nonlinear system of equations through Newton's method*)
m = 50;(*number of iterations*)
w = 0.005;(*relaxation parameter*)
init[0] = MapThread[{#1, #2} &, {Join[{ui}, Reverse[Table[((ui - uf)/(b - a)) (i - a) + uf, {i, a + h, b - h, h}]], {uf}], Join[{zi}, Reverse[Table[((zi - zf)/(b - a)) (i - a) + zf, {i, a + h, b - h, h}]], {zf}]}] // Flatten;(*Initial test values*)
For[j = 0, j <= m, j++, residuals = Table[{{resid1}, {resid2}} /. i -> j, {j, 2, n - 1}] // Flatten; 
DFxmat = sparse /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}; 
Residvec = Join[{0, 0}, residuals /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}, {residbound1, residbound2}/. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}]; 
init[j + 1] = init[j] + w LinearSolve[DFxmat, -Residvec]] // AbsoluteTiming

(*Residual tolerance of Dirichlet boundary at x=a PLUS interior points*)
ResidTol = Total[((Table[(Abs[resid1]) /. i -> j, {j, 2, n - 1}]) + (Table[(Abs[resid2]) /. i -> j, {j, 2, n - 1}])) /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}]/(2 n);
Print["Residual Tolerance = ", ResidTol]
Residual Tolerance = 0.00009350635229

(*Residual tolerance of Dirichlet boundary at x=a PLUS interior points PLUS Neumann boundary at x=b*)
ResidTol = Total[((Table[(Abs[resid1]) /. i -> j, {j, 2, n - 1}]) + (Table[(Abs[resid2]) /. i -> j, {j, 2, n - 1}]) + (Abs[residbound1]) + (Abs[residbound2])) /. {u[i_] :> init[j][[2 i - 1]], z[i_] :> init[j][[2 i]]}]/(2 n);
Print["Residual Tolerance = ", ResidTol]
Residual Tolerance = 0.02602675038

(*Gather data points*)
listu = Take[init[m], {1, -1, 2}];
listz = Take[init[m], {2, -1, 2}];
range = Range[a, b, h];
Listu = MapThread[{#1, #2} &, {range, listu}];
Listz = MapThread[{#1, #2} &, {range, listz}];
uapprox = Interpolation[Listu, InterpolationOrder -> 8, Method -> "Spline"];
zapprox = Interpolation[Listz, InterpolationOrder -> 8, Method -> "Spline"];

(*Plot*)
Plot[{uapprox[x], zapprox[x]}, {x, a, b}, Frame -> True, FrameLabel -> {"x"}, LabelStyle -> Directive[Black, 20], PlotStyle -> {{Green, Thickness[0.005]}, {Blue, Thickness[0.005]}}, PlotRange -> Full, ImageSize -> Large]

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