ClearAll["Global’*"] | FunM[fun_, X_]:= Module[{faux, dim, mataux, JordanD, sim, JordanF, eps, | fdiag, diagQ, fauxD}, (dim = Length@X; | faux[xx_, i_, j_]:= | Which[i <= j, 1/Abs[i - j]! (D[fun, {x, Abs[i - j]}]) /. x -> xx, True, 0]; | mataux[Y_]:= Table[faux[Y[[i, j]], i, j], {i, 1, dim}, {j, 1, dim}]; | JordanD = JordanDecomposition[X] // N; sim = JordanD[[]]; | JordanF = JordanD[[]]; eps = 1*10⋀-10; | diagQ = Norm[JordanF - DiagonalMatrix[Diagonal[JordanF]]]; | fauxD[xx_]:= (fun) /. x -> xx; | fdiag:= DiagonalMatrix[Map[fauxD, Diagonal[JordanF]]]; | Which[diagQ < eps, sim.fdiag.Inverse[sim], True, | sim.mataux[JordanF].Inverse[sim]])] | MGM[A_, B_, maxIterations_, tolerance_]:= Module[{k = 0}, | {n, n} = Dimensions[A]; Id = SparseArray[i_, i_} -> 1.}, {n, n}]; | A1 = SparseArray@ArrayFlatten[0, N@A}, {Id, 0]; Y[] = A1; | R[] = 1; Id2 = SparseArray[i_, i_} -> 1.}, {2 n, 2 n}]; | {Quiet@While[k < maxIterations && R[k] >= tolerance, | Y2 = Y[k].Y[k]; Y3 = Y2.Y[k]; Y4 = Y3.Y[k]; | Y5 = Y4.Y[k]; Y6 = Y5.Y[k]; Y7 = Y6.Y[k]; Y8 = Y7.Y[k]; | l1 = SparseArray[7 Id2 + 148 Y2 + 330 Y4 + 148 Y6 + 7 Y8]; | l2 = SparseArray@ArrayFlatten[Inverse@l1[[1;; n, 1;; n]], 0}, | {0, Inverse@l1[[n + 1;; 2 n, n + 1;; 2 n]]]; | Y[k + 1] = SparseArray[(48 Y[k] + 272 Y3 + 272 Y5 + 48 Y7).l2]; | R[k + 1] = Norm[Y[k + 1] - Y[k], Infinity]/ Norm[Y[k + 1], Infinity]; k++]; | AS = Y[k][[1;; n, n + 1;; 2 n]]; IAS = Y[k][[n + 1;; 2 n, 1;; n]]; | Mat = (IAS.B.IAS)}; | AS.FunM[Sqrt[x], Mat].AS] |
|