3
$\begingroup$

I had asked a question here regarding solution to a BVP problem. bbgodfrey provided an excellent answer using the method of integrated least squares.

However, for a few specific set of values of practical interest, the solution becomes unstable, especially at the boundaries.

Values for which the solution works perfectly:

L = 0.025; l = 0.025; w = 0.002; k = 390; A = 5625 10^-6; h = 3000; bh = 5.375; bc = 4 bh; λc = (1/l) (k w L/(A h/bc)); λh = (1/L) (k w l/(A h/bh)); V = 0.25;

Values for which the solution behaves wierdly:

L = 0.025; l = 0.025; w = 0.0001; k = 390; A = 5625 10^-6; h = 3000; bh = 5.375; bc = 4 bh; λc = (1/l) (k w L/(A h/bc)); λh = (1/L) (k w l/(A h/bh)); V = 0.25

Code

ClearAll[Evaluate[Context[] <> "*"]]
eq1 = D[θh[x, y], x] + bh (θh[x, y] - θw[x, y]);
eq2 = D[θc[x, y], y] + bc (θc[x, y] - θw[x, y]);
eq3 = λh D[θw[x, y], x, x] + λc V D[θw[x, y], y, y] + bh (θh[x, y] - θw[x, y]) + V bc (θc[x, y] - θw[x, y]);
th = Collect[(eq1 /. {θh -> Function[{x, y}, θhx[x] θhy[y]], θw -> Function[{x, y}, θwx[x] θwy[y]]})/(θhy[y] θwx[x]), {θhx[x], θhx'[x], θwy[y]}, Simplify];
1 == th[[1 ;; 3 ;; 2]];
eq1x = Subtract @@ Simplify[θwx[x] # & /@ %] == 0;
1 == -th[[2]];
eq1y = θhy[y] # & /@ %;
tc = Collect[(eq2 /. {θc -> Function[{x, y}, θcx[x] θcy[y]], θw -> Function[{x, y}, θwx[x] θwy[y]]})/(θcx[x] θwy[y]), {θcy[y], θcy'[y], θwy[y]}, Simplify];
1 == -tc[[1]];
eq2x = θcx[x] # & /@ %;
1 == tc[[2 ;; 3]];
eq2y = Subtract @@ Simplify[θwy[y] # & /@ %] == 0;
tw = Plus @@ ((List @@ Expand[eq3 /. {θh -> Function[{x, y}, θhx[x] θhy[y]], θc -> Function[{x, y}, θcx[x] θcy[y]], θw -> Function[{x, y}, θwx[x] θwy[y]]}])/(θwx[x] θwy[y]) /. Rule @@ eq2x /. Rule @@ eq1y);
sw == -tw[[1 ;; 5 ;; 2]];
eq3x = Subtract @@ Simplify[θwx[x] # & /@ %] == 0;
sw == tw[[2 ;; 6 ;; 2]];
eq3y = -Subtract @@ Simplify[θwy[y] # & /@ %] == 0;
sy = DSolveValue[{eq2y, eq3y, θcy[0] == 0, θwy'[0] == 0}, {θwy[y], θcy[y], θwy'[1]}, {y, 0, 1}] /. C[2] -> 1;
sx = DSolveValue[{eq1x, eq3x, θwx'[0] == 0, θwx'[1] == 0, θhx[0] == 1}, {θwx[x], θhx[x]}, {x, 0, 1}];

L = 0.025; l = 0.025; w = 0.0001; k = 390; A = 5625 10^-6; h = 3000; bh = 5.375; bc = 4 bh; λc = (1/l) (k w L/(A h/bc)); λh = (1/L) (k w l/(A h/bh)); V = 0.25;
disp = sy[[3]];
n = 12;
Plot[disp, {sw, -800, 10}, AxesLabel -> {sw, "disp"}, LabelStyle -> {15, Bold, Black}, ImageSize -> Large]
Partition[Union @@ Cases[%, Line[z_] -> z, Infinity], 2, 1];
Reverse[Cases[%, {{z1_, z3_}, {z2_, z4_}} /; z3 z4 < 0 :> z1]][[1 ;; n]];
tsw = sw /. Table[FindRoot[disp, {sw, sw0}], {sw0, %}];
syn = ComplexExpand@Replace[bh sy[[1]] /. C[2] -> 1, {sw -> #} & /@ tsw, Infinity] //Chop // Chop;
Integrate[Expand[(1 - Array[c, n].syn)^2], {y, 0, 1}] // Chop // Chop;
coef = NArgMin[%, Array[c, n]]
Plot[coef.syn - 1, {y, 0, 1}, AxesLabel -> {y, err}, LabelStyle -> {15, Bold, Black}, PlotRange -> Full, ImageSize -> Large]
solw = coef.ComplexExpand@Replace[sy[[1]] sx[[1]], {sw -> #} & /@ tsw, Infinity];
solh = coef.ComplexExpand@Replace[bh sy[[1]] sx[[2]], {sw -> #} & /@ tsw, Infinity];
solc = coef.ComplexExpand@Replace[bc sy[[2]] sx[[1]], {sw -> #} & /@ tsw, Infinity];
cavg = Integrate[solc, {x, 0, 1}] // Chop;
havg = Integrate[solh, {y, 0, 1}] // Chop;
{cavg /. y -> 1, havg /. x -> 1}

{Plot3D[solc, {x, 0, 1}, {y, 0, 1}], Plot3D[solh, {x, 0, 1}, {y, 0, 1}]}

It is pretty much evident that the problem occurs due to very closely spaced roots (eigen-values) for the second set of constants. However, I am not sure if I am catching it correctly.

Improve this question
$\endgroup$
4
  • $\begingroup$ Do you try to solve PVP or to make research about stability of least squares algorithm? $\endgroup$
    – Alex Trounev
    Jan 21, 2022 at 14:39
  • $\begingroup$ @AlexTrounev My objective was to solve this BVP using an anlytical/semi-analytical approach. The linked question has the original answer which proposed integrated least squares as an alternative to do this. The parameter sets I have written here are just describing two different flow configurations. Stability of least squares algorith is not my area of research or interest. $\endgroup$
    – Avrana
    Jan 21, 2022 at 15:03
  • $\begingroup$ I have solution for first and second set of parameters based on Euler wavelets. For the first set the difference between my solution and yours is about 10^-4. In a case of second set of parameters solution looks same as for the first case. So, I don't understand why we can't compute this solution with the least squares algorithm. $\endgroup$
    – Alex Trounev
    Jan 23, 2022 at 17:23
  • $\begingroup$ @AlexTrounev Thankyou for the effort Alex. However, I am afraid I am unfamiliar with the method of Euler wavelets. Please post your solution if possible. Also, were you able to get a well-behaving solution at the boundaries with my code ? Additionally, the parameter k=390 represents copper (thermal conductivity). However, when I switch to steel, i.e. k=16, my method always behaves unbounded at the boundaries. $\endgroup$
    – Avrana
    Jan 23, 2022 at 17:44

1 Answer 1

Reset to default
7
+50
$\begingroup$

We can solve this problem with using numerical method developed in our recent paper and based on Euler wavelets. First we define wavelets and all functions to be computed

UE[m_, t_] := EulerE[m, t]
psi[k_, n_, m_, t_] := 
 Piecewise[{{2^(k/2) Sqrt[2/Pi] UE[m, 2^k t - 2 n + 1], (n - 1)/
      2^(k - 1) <= t < n/2^(k - 1)}, {0, True}}]
PsiE[k_, M_, t_] := 
 Flatten[Table[psi[k, n, m, t], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]
k0 = 3; M0 = 4; With[{k = k0, M = M0}, 
 var1 = Flatten[Table[c[n, m], {n, 1, 2^(k - 1)}, {m, 0, M - 1}]]];
nn = Length[var1];
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; ycol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 With[{k = k0, M = M0}, PsiE[k, M, t1]]; Int1 = 
 With[{k = k0, M = M0}, Integrate[PsiE[k, M, t1], t1]];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y; 
int2[y_] := Int2 /. t1 -> y; M = nn;

U1 = Array[a1, {M, M}]; U2 = Array[a2, {M, M}]; U3 = 
 Array[a3, {M, M}]; U4 = Array[a4, {M, M}]; G1 = Array[g1, {M}]; G2 = 
 Array[g2, {M}]; G3 = Array[g3, {M}]; G4 = Array[g4, {M}]; F1 = 
 Array[f1, {M}]; F2 = Array[f2, {M}];


u1[x_, y_] := int2[x] . U1 . Psi[y] + x G1 . Psi[y] + F1 . Psi[y]; 
u2[x_, y_] := Psi[x] . U2 . int2[y] + y G2 . Psi[x] + F2 . Psi[x];
uy[x_, y_] := Psi[x] . U2 . int1[y] + G2 . Psi[x];
ux[x_, y_] := int1[x] . U1 . Psi[y] + G1 . Psi[y];
uxx[x_, y_] := Psi[x] . U1 . Psi[y];
uyy[x_, y_] := Psi[x] . U2 . Psi[y];
tehx[x_, y_] := Psi[x] . U3 . Psi[y]; 
tecy[x_, y_] := Psi[x] . U4 . Psi[y]; 
teh[x_, y_] := int1[x] . U3 . Psi[y] + G3 . Psi[y]; 
tec[x_, y_] := Psi[x] . U4 . int1[y] + G4 . Psi[x];

Please, note, that we solve system of equations eq1 = D[θh[x, y], x] + bh (θh[x, y] - θw[x, y]); eq2 = D[θc[x, y], y] + bc (θc[x, y] - θw[x, y]); eq3 = λh D[θw[x, y], x, x] + λc V D[θw[x, y], y, y] + bh (θh[x, y] - θw[x, y]) + V bc (θc[x, y] - θw[x, y]);, and u1, teh, tec are θw, θh, θc consequently. Second step, we define system of equations and variables on the grid for k=16, w=0.0001. Test 1:

L = 0.025; l = 0.025; w = 0.0001; k = 16; A = 
 5625 10^-6; h = 3000; bh = 5.375; bc = 
 4 bh; \[Lambda]c = (1/l) (k w L/(A h/bc)); \[Lambda]h = (1/
    L) (k w l/(A h/bh)); V = 0.25; var = 
 Join[Flatten[U1], Flatten[U2], G1, G2, F1, F2, Flatten[U3], 
  Flatten[U4], G3, G4]; eq = 
 Join[Flatten[
   Table[\[Lambda]h*uxx[xcol[[i]], ycol[[j]]] + \[Lambda]c*V*
       uyy[xcol[[i]], ycol[[j]]] - tehx[xcol[[i]], ycol[[j]]] - 
      V*tecy[xcol[[i]], ycol[[j]]] == 0, {i, M}, {j, M}]], 
  Flatten[Table[
    tehx[xcol[[i]], ycol[[j]]] + 
      bh (teh[xcol[[i]], ycol[[j]]] - u1[xcol[[i]], ycol[[j]]]) == 
     0, {i, M}, {j, M}]], 
  Flatten[Table[
    tecy[xcol[[i]], ycol[[j]]] + 
      bc (tec[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]]) == 
     0, {i, M}, {j, M}]], 
  Flatten[Table[
    u1[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]] == 0, {i, 
     M}, {j, M}]], 
  Flatten[Table[{ux[1, ycol[[j]]] == 0, ux[0, ycol[[j]]] == 0, 
     uy[xcol[[j]], 0] == 0, uy[xcol[[j]], 1] == 0, 
     teh[0, ycol[[j]]] == 1, tec[xcol[[j]], 0] == 0}, {j, M}]]];

Finally we calculate and visualize numerical solution as follows

sol = FindRoot[eq, Table[{var[[i]], 1/10}, {i, Length[var]}], 
   MaxIterations -> 1000];

{Plot3D[Evaluate[teh[x, y] /. sol], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]h, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White], 
 Plot3D[Evaluate[tec[x, y] /. sol], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]c, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White], 
 Plot3D[Evaluate[u1[x, y] /. sol], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]w, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White]}

Figure 1

Test 2. We compute solution for w = 0.0001; k = 390; Figure 2

Test 3. We compute solution for w = 0.002; k = 390;. In this case we can compare solution computed with wavelets - Figure 3 and with least squares algorithm - Figure 4. In the last picture shown difference two solutions θw on the grid Figure 3 Figure 4

Update 1. Implementation of list squares method as it explained in the paper LEAST SQUARES FOURIER SERIES SOLUTIONS TO BOUNDARY VALUE PROBLEMS by ROBERT B. KELMAN. Definitions of base function and solution

UE[m_, t_] := Cos[m t] Exp[-m t]
nn = 5;
dx = 1/(nn);  xl = Table[ l*dx, {l, 0, nn}]; ycol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 Table[UE[n , t1], {n, 0, nn - 1}]; Int1 = Integrate[Psijk, t1];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y; 
int2[y_] := Int2 /. t1 -> y; M = nn;

U1 = Array[a1, {M, M}]; U2 = Array[a2, {M, M}]; U3 = 
 Array[a3, {M, M}]; U4 = Array[a4, {M, M}]; G1 = Array[g1, {M}]; G2 = 
 Array[g2, {M}]; G3 = Array[g3, {M}]; G4 = Array[g4, {M}]; F1 = 
 Array[f1, {M}]; F2 = Array[f2, {M}];


u1[x_, y_] := int2[x] . U1 . Psi[y] + x G1 . Psi[y] + F1 . Psi[y]; 
u2[x_, y_] := Psi[x] . U2 . int2[y] + y G2 . Psi[x] + F2 . Psi[x];
uy[x_, y_] := Psi[x] . U2 . int1[y] + G2 . Psi[x];
ux[x_, y_] := int1[x] . U1 . Psi[y] + G1 . Psi[y];
uxx[x_, y_] := Psi[x] . U1 . Psi[y];
uyy[x_, y_] := Psi[x] . U2 . Psi[y];
tehx[x_, y_] := Psi[x] . U3 . Psi[y]; 
tecy[x_, y_] := Psi[x] . U4 . Psi[y]; 
teh[x_, y_] := int1[x] . U3 . Psi[y] + G3 . Psi[y]; 
tec[x_, y_] := Psi[x] . U4 . int1[y] + G4 . Psi[x];

We solve system of equations eq1 = D[θh[x, y], x] + bh (θh[x, y] - θw[x, y]); eq2 = D[θc[x, y], y] + bc (θc[x, y] - θw[x, y]); eq3 = λh D[θw[x, y], x, x] + λc V D[θw[x, y], y, y] + bh (θh[x, y] - θw[x, y]) + V bc (θc[x, y] - θw[x, y]);, and u1, teh, tec are θw, θh, θc consequently. Definitions system of equations and variables on the grid for k=16, w=0.0001. Test 1:

L = 0.025; l = 0.025; w = 0.0001; k = 16; A = 
 5625 10^-6; h = 3000; bh = 5.375; bc = 
 4 bh; \[Lambda]c = (1/l) (k w L/(A h/bc)); \[Lambda]h = (1/
    L) (k w l/(A h/bh)); V = 0.25; var = 
 Join[Flatten[U1], Flatten[U2], G1, G2, F1, F2, Flatten[U3], 
  Flatten[U4], G3, G4]; eq = 
 Join[Flatten[
   Table[\[Lambda]h*uxx[xcol[[i]], ycol[[j]]] + \[Lambda]c*V*
       uyy[xcol[[i]], ycol[[j]]] - tehx[xcol[[i]], ycol[[j]]] - 
      V*tecy[xcol[[i]], ycol[[j]]] == 0, {i, M}, {j, M}]], 
  Flatten[Table[
    tehx[xcol[[i]], ycol[[j]]] + 
      bh (teh[xcol[[i]], ycol[[j]]] - u1[xcol[[i]], ycol[[j]]]) == 
     0, {i, M}, {j, M}]], 
  Flatten[Table[
    tecy[xcol[[i]], ycol[[j]]] + 
      bc (tec[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]]) == 
     0, {i, M}, {j, M}]], 
  Flatten[Table[
    u1[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]] == 0, {i, 
     M}, {j, M}]], 
  Flatten[Table[{ux[1, ycol[[j]]] == 0, ux[0, ycol[[j]]] == 0, 
     uy[xcol[[j]], 0] == 0, uy[xcol[[j]], 1] == 0, 
     teh[0, ycol[[j]]] == 1, tec[xcol[[j]], 0] == 0}, {j, M}]]];

Finally we calculate and visualize numerical solution

{bvec, mat} = CoefficientArrays[eq, var];

sol = LinearSolve[mat, -bvec];

rul = Table[var[[i]] -> sol[[i]], {i, Length[var]}];

{Plot3D[Evaluate[teh[x, y] /. rul], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]h, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White], 
 Plot3D[Evaluate[tec[x, y] /. rul], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]c, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White], 
 Plot3D[Evaluate[u1[x, y] /. rul], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]w, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White]}

Figure 5

Update 2. We always can solve this problem as an optimization (minimization) problem as follows

UE[m_, t_] := Cos[m t] Exp[-m t]
nn = 5;
dx = 1/(nn); xl = Table[l*dx, {l, 0, nn}]; ycol = 
 xcol = Table[(xl[[l - 1]] + xl[[l]])/2, {l, 2, nn + 1}]; Psijk = 
 Table[UE[n, t1], {n, 0, nn - 1}]; Int1 = Integrate[Psijk, t1];
Int2 = Integrate[Int1, t1];
Psi[y_] := Psijk /. t1 -> y; int1[y_] := Int1 /. t1 -> y;
int2[y_] := Int2 /. t1 -> y; M = nn;

U1 = Array[a1, {M, M}]; U2 = Array[a2, {M, M}]; U3 = 
 Array[a3, {M, M}]; U4 = Array[a4, {M, M}]; G1 = Array[g1, {M}]; G2 = 
 Array[g2, {M}]; G3 = Array[g3, {M}]; G4 = Array[g4, {M}]; F1 = 
 Array[f1, {M}]; F2 = Array[f2, {M}];


u1[x_, y_] := int2[x] . U1 . Psi[y] + x G1 . Psi[y] + F1 . Psi[y];
u2[x_, y_] := Psi[x] . U2 . int2[y] + y G2 . Psi[x] + F2 . Psi[x];
uy[x_, y_] := Psi[x] . U2 . int1[y] + G2 . Psi[x];
ux[x_, y_] := int1[x] . U1 . Psi[y] + G1 . Psi[y];
uxx[x_, y_] := Psi[x] . U1 . Psi[y];
uyy[x_, y_] := Psi[x] . U2 . Psi[y];
tehx[x_, y_] := Psi[x] . U3 . Psi[y];
tecy[x_, y_] := Psi[x] . U4 . Psi[y];
teh[x_, y_] := int1[x] . U3 . Psi[y] + G3 . Psi[y];
tec[x_, y_] := Psi[x] . U4 . int1[y] + G4 . Psi[x];
L = 0.025; l = 0.025; w = 0.0001; k = 16; A = 
 5625 10^-6; h = 3000; bh = 5.375; bc = 
 4 bh; \[Lambda]c = (1/l) (k w L/(A h/bc)); \[Lambda]h = (1/
    L) (k w l/(A h/bh)); V = 0.25; var = 
 Join[Flatten[U1], Flatten[U2], G1, G2, F1, F2, Flatten[U3], 
  Flatten[U4], G3, G4]; eq = 
 Join[Flatten[
   Table[\[Lambda]h*uxx[xcol[[i]], ycol[[j]]] + \[Lambda]c*V*
      uyy[xcol[[i]], ycol[[j]]] - tehx[xcol[[i]], ycol[[j]]] - 
     V*tecy[xcol[[i]], ycol[[j]]], {i, M}, {j, M}]], 
  Flatten[Table[
    tehx[xcol[[i]], ycol[[j]]] + 
     bh (teh[xcol[[i]], ycol[[j]]] - u1[xcol[[i]], ycol[[j]]]), {i, 
     M}, {j, M}]], 
  Flatten[Table[
    tecy[xcol[[i]], ycol[[j]]] + 
     bc (tec[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]]), {i, 
     M}, {j, M}]], 
  Flatten[Table[
    u1[xcol[[i]], ycol[[j]]] - u2[xcol[[i]], ycol[[j]]], {i, M}, {j, 
     M}]]]; bc = 
 Flatten[Table[{ux[1, ycol[[j]]] == 0, ux[0, ycol[[j]]] == 0, 
    uy[xcol[[j]], 0] == 0, uy[xcol[[j]], 1] == 0, 
    teh[0, ycol[[j]]] == 1, tec[xcol[[j]], 0] == 0}, {j, M}]];

sol = NMinimize[{eq . eq, bc}, var]

Visualization

{Plot3D[Evaluate[teh[x, y] /. sol[[2]]], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]h, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White], 
 Plot3D[Evaluate[tec[x, y] /. sol[[2]]], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]c, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White], 
 Plot3D[Evaluate[u1[x, y] /. sol[[2]]], {x, 0, 1}, {y, 0, 1}, 
  PlotRange -> All, ColorFunction -> "Rainbow", 
  AxesLabel -> Automatic, PlotLabel -> \[Theta]w, PlotPoints -> 50, 
  PlotTheme -> "Scientific", MeshStyle -> White]}

Figure 6

Improve this answer
$\endgroup$
16
  • $\begingroup$ Thankyou for the answer. I will initially go through the attached paper to get some idea about the method and come back to this answer, to properly understand it. $\endgroup$
    – Avrana
    Jan 24, 2022 at 5:22
  • $\begingroup$ The referred paper is excellent (absolute errors are basicallly null with the method). Although, I still need to read more to properly understand this approach. However, I have checked that the code works for all material (k=16, k=390) and operating conditions (different bh,bc, lambdac, lambdah, V). Hence, I will mark this answer as accepted. Although, I still do not understand why the least squares code could not converge properly for k=16 configurations. Any questions, I have regarding your paper, hopefully, I can enquire along this thread. $\endgroup$
    – Avrana
    Jan 25, 2022 at 5:24
  • 1
    $\begingroup$ @Avrana We can extend this method for 3D case as well that is not discussed in the paper linked. Also we can improve code proposed by bbgodfrey to make it more stable in desired parameters range. $\endgroup$
    – Alex Trounev
    Jan 25, 2022 at 7:06
  • 1
    $\begingroup$ @Avrana We always can solve this problem as minimization problem, but there is no fast solver for numerical minimization, and it is why we use LinearSolve. See Update 2 with NMinimize as solver. Equation u1=u2 on the grid means that $\int(\int u_{xx}dx)dx=\int(\int u_{yy}dy)dy$. Yes, I have my code with wavelets as a step to develop last code. $\endgroup$
    – Alex Trounev
    Jan 27, 2022 at 10:54
  • 1
    $\begingroup$ @Avrana See my post on community.wolfram.com/groups/-/m/t/2457908 $\endgroup$
    – Alex Trounev
    Jan 28, 2022 at 6:52

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.