# Updating a parameter each time step within a Finite Difference scheme?

So I am solving a PDE for a function $h(\theta,t)$ via finite difference scheme. The PDE has a function $Q$ in it, which I wish to update each time step depending on where $h$ lies.

Firstly, let us look at simple case when $Q = const$. Here is my code:

ClearAll[dh, t, n, ds, h0, sol, h, valh, ti, s];

n = 1001;
ds = Pi/(n - 1);

(* dh: R^n \[Rule] R^n, symbolic and numeric *)
dh[h_List /; Length@h == n] := With[{

s = Array[# &, n, {-π, 0}],
Q = 0.005,(*I wish to update this each time step depending on \
where h lies.*)
d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #] &,
d2 = ListCorrelate[{1, -2, 1}/ds^2, #] &,
h2 = Join[{2}, h, {2}]
}, (-2 + h) (-1 +
h) (-18 Q d1[h2] + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[
s] d1[h2] + (-2 + h) (-1 +
h) ((-9 + 6 h) d1[h2]^2 - (2 - 3 h + h^2) (h - d2[h2])) Sin[
s])
];

(*h0 ∈ R^n *)
h0 = Table[Piecewise[{

{1.999, -Pi < θ < -7 Pi/8},
{(2 Sin[-6 Pi/7])/
Sin[θ], -6 Pi/
7 < θ < -2.0909003694838066}, {1.001, \
-2.26998149133592 < θ < -0.6662394324925154},
{(2 Sin[-(Pi/10)])/
Sin[θ], -0.6662394324925154 < θ < -(Pi/10)},
{1.999, -(Pi/8) < θ < 0}

}, 1.999]

, {θ, -1. Pi, 0., ds}
];

(*Solve
dh/dt(t) = dh(t)
h(0)= h0

for h(t) ∈ R^n
for t ∈ (0,20)
*)
sol = With[{
h = h[#][t] & /@ Range@n
},

NDSolveValue[
{
D[h, t] == dh[h],
(h /. t -> 0) == h0
},
{t, 0, 20}
]

];

valh = Transpose[#["ValuesOnGrid"] & /@ sol];
ti = Flatten@sol[[1]]["Grid"];
s = Array[# &, n, {-Pi, 0}];

h = ListInterpolation[valh, {ti, s}, InterpolationOrder -> 1];

ListAnimate@Table[
Show[
(*Interface*)
PolarPlot[
Evaluate[h[ti, θ]],
{θ, -π, 0},

PlotRange -> {{-2, 2}, {-2, 2}},
ImageSize -> Large,
PlotStyle -> Black],
(*Ring*)
RegionPlot[
1 <= x^2 + y^2 <= 4,
{x, -2, 2}, {y, -2, 2},

PlotStyle -> Opacity[0.1], AspectRatio -> 1
]

]
, {ti, 0, 20, 0.2}]


As we can see the interface moves around with a fixed motion, due to Q being fixed at 0.005. However, in reality I wish to update $Q$ at each time step with the following relation $Q = \cos \psi - \cos \phi$, where $\phi$ and $\psi$ are the angles from the base of the annulus to the middle of the right and left flat sections.

Now I use the following pattern matching code to find the initial $Q$. It basically looks for the median position of the centre of the interface and multiply by ds to convert it to an angle.

Q = N@Cos[Median@Position[h0[[1 ;; 500]], x_ /; 1.001 < x < 1.999] ds][[1]] - N@Cos[(Median@
Position[h0[[501 ;; 1001]], x_ /; 1.001 < x < 1.999]) ds][[1]]


This gives me a sensible initial $Q$ of 0.26. I wish to then update this $Q$ with each time step, with the new position of the middle of the interface. My hope is eventually that the two angles will become equal and $Q = 0$. Here is my attempt to update $Q$ at each time step in the Finite Difference scheme.

n = 1001;
ds = Pi/(n - 1);
ClearAll[dh, t];
dh[h_List] :=
With[{s = Array[# &, n, {-π, 0}],
Q = N@Cos[Median@Position[#[[1 ;; 500]], x_ /; 1.001 < x < 1.999] ds][[1]]-N@Cos[(Median@Position[#[[501 ;; 1001]], x_ /; 1.001 < x < 1.999]) ds[[1]] &,
d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #] &,
d2 = ListCorrelate[{1, -2, 1}/ds^2, #] &},
(-2 + h) (-1 + h) (-18 Q[#1] d1[#1] + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[
s] d1[#1] + (-2 + h) (-1 +
h) ((-9 + 6 h) d1[#1]^2 - (2 - 3 h + h^2) (h -
d2[#1])) Sin[s]) &@Join[{2}, h, {2}]];


Which runs, but nothing moves! It should really as the initial $Q$ for $h_0$ is 0.26 as shown before.

Any pointers on how to update $Q$ at each time step? I feel like I'm nearly there but perhaps my syntax is incorrect? Maybe using the pure function incorrectly?

Any help appreciated as always.

Edit

So another idea I had was to include $Q$ as an unknown in the finite difference scheme. ie. It would appear as Q[1], Q[2] etc. in each of the equations. I could then solve the finite difference scheme D[h,t] == dh[h,Q] coupled with an equation for $Q$, in the NDSolveValue. Here is my attempt:

n = 1001;
ds = Pi/(n - 1);
ClearAll[dh, t];
dh[h_List, Q_] :=
With[{s = Array[# &, n, {-π, 0}],
d1 = ListCorrelate[{-0.5, 0, 0.5}/ds, #] &,
d2 = ListCorrelate[{1, -2, 1}/ds^2, #] &},
(-2 + h) (-1 +
h) (-18 Q d1[#1] + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[
s] d1[#1] + (-2 + h) (-1 +
h) ((-9 + 6 h) d1[#1]^2 - (2 - 3 h + h^2) (h -
d2[#1])) Sin[s]) &@Join[{2}, h, {2}]];

q0 = 0.26

Clear@h;
sol = With[{h = h[#][t] & /@ Range@n, Q = Q[#][t] & /@ Range@n},
NDSolveValue[{D[h, t] == dh[h, Q],
Q[t] == -N@
Cos[(Median@
Position[h[[501 ;; 1001]],
x_ /; 1.001 < x < 1.999]) ds][[1]] +
N@Cos[Median@
Position[h[[1 ;; 500]],
x_ /; 1.001 < x < 1.999] ds][[1]] , (h /. t -> 0) ==
h0, (Q /. t -> 0) == q0}, Head /@ h, {t, 0, 20},
Method -> {"EquationSimplification" -> "Residual"}]];


But this gives me a very strange error.

Transpose::nmtx: The first two levels of {NDSolvexs$422676, NDSolvexs$422677, NDSolvexs$422678, NDSolvexs$422679, NDSolvexs$422680, NDSolvexs$422681,NDSolvexs$422682, NDSolvexs$422683, NDSolvexs$422684, etc. <<951>>} cannot be transposed. • Perhaps, the option StepMonitor would be helpful. – bbgodfrey Aug 3 '16 at 14:14 • @bbgodfrey Does look like it could be useful option. Not quite sure how to implement it into the NDSolveValue? – M Horsley Aug 3 '16 at 15:46 • I have not used StepMonitor, so I cannot offer practical advice. If you search Mathematica.Stackexchange and also the web with Google, you will find several examples and a few references. When you solve this problem, please be sure to provide it as an answer to your question. Many others may benefit. – bbgodfrey Aug 3 '16 at 15:58 • @bbgodfrey Using your comment as inspiration, if I add a Sow tag to the finite difference scheme, and then a Reap tag around NDSolveValue. Then sol[[2,1,1,1]] returns the first equation from the scheme. It seems that$Q$is 0 in this case, and continues to be. This would explain why nothing moves. So there's clearly a problem with my pure function syntax for Q as it should be passing 0.26 in the first equation. The syntax works for$h_0$case so I believe the pattern matching bit is fine. It's just how to pass each successive$h$into the$Q$calculation. – M Horsley Aug 3 '16 at 16:30 • "where$ϕ$and$ψ$are the angles from the base of the annulus to the middle of the right and left flat sections", then I think your definition for Q in the code is wrong, notice e.g. h0[[1 ;; 500]] isn't just formed by points on the flat section. Also, as shown in the .gif, when the section begins to move, it's no longer flat, how will you decide the middle in this case? – xzczd Aug 5 '16 at 3:21 ## 2 Answers OK, finally find some time to write the answer. First I'd like to point out why your 2 approaches failed. As mentioned in the comment above, they failed for the same reason: NDSolve doesn't have Hold* attribute, the Position in your code evaluates before the equation is sent into NDSolve. For more information you may have a look at this post about essentially the same problem. Then I'd like to point out that your definition for$\psi$in the code is wrong, there should be a Pi/2 - in the first Cos[…]. Finally the following is my solution. It's based on Michael E2's answer in your previous question: ClearAll[dh2, s, h]; n = 201; s = N@Array[# &, n, {- π, 0}]; ds = Pi/(n - 1); dhC = With[{s = s, findmid = Last@Median@Select[#, 1.001 < #[[1]] < 1.999 &] &}, Compile[{{h, _Real, 1}, {d1, _Real, 1}, {d2, _Real, 1}}, Module[{halfl = Floor[Length@h/2], Q}, Q = Cos[Pi/2 + findmid[{h, s}\[Transpose][[;; halfl]]]] - Cos[-Pi/2 - findmid[{h, s}\[Transpose][[halfl + 1 ;; -1]]]]; (-2 + h) (-1 + h) (-18 Q d1 + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[ s] d1 + (-2 + h) (-1 + h) ((-9 + 6 h) d1^2 - (2 - 3 h + h^2) (h - d2)) Sin[ s])](*, CompilationTarget -> C*)]]; With[{(*compute derivatives and Sin[s],Cos[s] just once*) d1FN = NDSolveFiniteDifferenceDerivative[1, s, "DifferenceOrder" -> 2], d2FN = NDSolveFiniteDifferenceDerivative[2, s, "DifferenceOrder" -> 1]}, dh2[h_List] := With[{d1 = d1FN@h, d2 = d2FN@h}, dhC[h, d1, d2]]]; h0 = Table[Piecewise[{{1.999, -Pi < θ < -7 Pi/8}, {(2 Sin[-6 Pi/7])/Sin[θ], -6 Pi/7 < θ < -2.0909003694838066}, {1.001, -2.26998149133592 < θ < -0.6662394324925154}, {(2 Sin[-(Pi/10)])/Sin[θ], -0.6662394324925154 < θ < -(Pi/10)}, {1.999, -(Pi/8) < θ < 0}}, 1.999], {θ, -1. Pi, 0., ds}]; {sol2} = h /. NDSolve[{h'[t] == dh2[h[t]], h[0] == h0}, h, {t, 0, 8.43}, MaxSteps -> Infinity]; // AbsoluteTiming valh = sol2["ValuesOnGrid"]; ti = Flatten@sol2["Grid"]; newsol = ListInterpolation[valh, {ti, s}] pic = RegionPlot[1 <= x^2 + y^2 <= 4, {x, -2, 2}, {y, -2, 2}, PlotStyle -> Opacity[0.1]]; (*ListAnimate@*)Export["a.gif", Table[Show[PolarPlot[newsol[t, θ], {θ, -Pi, 0}, PlotRange -> 2, PlotStyle -> Black], pic], {t, 0, 8.4, 0.1}]]  I modified the end time to 8.43 because when the calculation got stuck at this point, I'm not sure about the reason, but relatively confident about the coding, perhaps it's the nature of the model? Update OP asked if it's possible to monitor the variation of Q, the answer is of course yes. We just need to add a Sow@Q in dhC: dhC = With[{s = s, findmid = Last@Median@Select[#, 1.001 < #[[1]] < 1.999 &] &}, Compile[{{h, _Real, 1}, {d1, _Real, 1}, {d2, _Real, 1}}, Module[{halfl = Floor[Length@h/2], Q}, Q = Cos[Pi/2 + findmid[{h, s}\[Transpose][[;; halfl]]]] - Cos[-Pi/2 - findmid[{h, s}\[Transpose][[halfl + 1 ;; -1]]]]; Sow@Q; (-2 + h) (-1 + h) (-18 Q d1 + (-2 + h) (-1 + h) (4 + h (-15 + 8 h)) Cos[ s] d1 + (-2 + h) (-1 + h) ((-9 + 6 h) d1^2 - (2 - 3 h + h^2) (h - d2)) Sin[ s])](*,CompilationTarget\[Rule]C*)]];  Notice Sow is actually not compilable, but usually it won't slow down CompiledFunction that much. Here's another example. Then add a Reap outside of NDSolve: {{sol2}, {Qlst}} = Reap[h /. NDSolve[{h'[t] == dh2[h[t]], h[0] == h0}, h, {t, 0, 8.43}, MaxSteps -> Infinity]]; // AbsoluteTiming  OK, let's check the variation of Q: ListPlot[Qlst, PlotRange -> All]  It goes to 0, as you expected. • Hmm, I don't understand why the interface moves so much. The$Q$should tend to 0 and then stay there. Almost like John Bishop's solution does. You can see how the$Q$varies and goes to 0 - his goes beyond 0 which is another sauce of confusion for me. Is it possible to list the$Q$at each time step in your solution? StepMontior I guess in NDSolveValue ? – M Horsley Aug 10 '16 at 19:59 • @ojlm As mentioned in the comment under Bishop's answer, I think his WhenEvent doesn't work. As to the monitoring of Q, the answer is of course yes. See my edit. – xzczd Aug 11 '16 at 2:40 • Thanks for your update. It's great that we can see that$Q$goes to 0, but it goes to 0 when clearly the angles are not equal? Given the definition,$Q = 0$is when$\sin \psi = \sin \phi$. However, it's clear to me from the animation that when$Q = 0$we have$\phi \neq \psi$. – M Horsley Aug 11 '16 at 8:36 • @ojlm It's not surprising to me, because your algorithm for selecting$ψ$and$ϕ$is rather coarse, it takes the Median of points between 1.001<h<1.999, which involves a lot of undesired points. (See those points which are quite close to the wall but clearly in the range of 1.001<h<1.999?) I believe using a better standard e.g. 1.4<h<1.6 will result in a better result, but with this standard the calculation will become so memory-consuming that my laptop can't bear it for a long enough time... – xzczd Aug 11 '16 at 9:42 I think that the meta-problem here is that MMA doesn't allow vector-valued functions in NDSolve (and related solvers). Given this restriction, I couldn't figure out how to define a real-valued function$Q(h(\theta,t))$that would incorporate into NDSolve with a vector set of equations$h(\theta_i,t)$, where$\theta_i\in\{x\times\pi/1000|{-1000}\leq x \leq 0\wedge x\in\mathbb{Z}\}$. (Does that sound right? I think it does.) Instead, I used WhenEvent inside NDSolveValue to discretely update$Q$as$Q(\theta_i,h(\theta_i,t))$for the set of equations$h(\theta_i,t)$at each time point$t=t_i\in\{x\times0.2|0\leq x \leq 100\wedge x\in\mathbb{Z}\}$. n = 1001; dθ = π/(n - 1); θ = Table[θ, {θ, -π, 0, dθ}]; functionList[h_, θ_List, t_: 0] := Table[h[s][t], {s, θ}]; h0 = N@Table[ Piecewise[{ {1.999, -π < θ < -((7 π)/8)}, {(2 Sin[-((6 π)/7)])/ Sin[θ], -((6 π)/ 7) < θ < -2.0909003694838066}, {1.001, -2.26998149133592 < θ < -0.6662394324925154}, {(2 Sin[-(π/10)])/ Sin[θ], -0.6662394324925154 < θ < -(π/10)}, {1.999, -(π/8) < θ < 0} }, 1.999], {θ, θ}]; dh[h_, θ_List, t_: 0, Q_: 0.005] := Module[{hθt = functionList[h, θ, t], d1, d2}, d1 = ListCorrelate[{-0.5, 0, 0.5}/dθ, Join[{2}, hθt, {2}]]; d2 = ListCorrelate[{1, -2, 1}/dθ^2, Join[{2}, hθt, {2}]]; (-2 + hθt) (-1 + hθt) (-18 Q d1 + (-2 + hθt) (-1 + hθt) (4 + hθt (-15 + 8 hθt)) Cos[θ] d1 + (-2 + hθt) (-1 + hθt) ((-9 + 6 hθt) d1^2 - (2 - 3 hθt + hθt^2) (hθt - d2)) Sin[θ]) ]  I defined the center-finding function as in the OP, but to return polar coordinates instead of just angles. Then I defined a helper function to find the interface center points for a given$h(\theta_i,t_i)$. Finally I defined$Q(\theta_i,h(\theta_i,t_i))$and a polymorphic shortcut to work on interface data lists. interfaceCenters[θ_List, hθt_List] := { {π/2 -dθ First@Median[#], Extract[hθt, Round[Median[#], 1]]} &@ Position[hθt[[1 ;; 500]], x_ /; 1.001 < x < 1.999], {dθ First@Median[#], Extract[hθt, Round[Median[#], 1] + 500]} &@ Position[hθt[[501 ;; 1001]], x_ /; 1.001 < x < 1.999] }; Q[θ_List, hθt_List] := With[ {iCs = interfaceCenters[θ, hθt]}, N@Cos[iCs[[1, 1]]] - N@Cos[iCs[[2, 1]]] ]; Q[iC_List] := N@Cos[#[[1, 1]]] - Cos[#[[2, 1]]] &@iC; Q0 = Q[\[Theta], h0];  Then the simulation and data generation: Clear[h]; sol = NDSolveValue[{ D[functionList[h, θ, t], t] == dh[h, θ, t,Q0], functionList[h, θ, 0] == h0, WhenEvent[Mod[t, 0.2] == 0, With[ {q = Q[θ, functionList[h, θ, t]]}, dh[h, θ, t,_Real] -> dh[h, θ, t, q]]; ] }, Head /@ functionList[h, θ, t], {t, 0, 20}]; valh = Transpose[#["ValuesOnGrid"] & /@ sol]; ti = Flatten@sol[[1]]["Grid"]; hSol = ListInterpolation[valh, {ti, θ}, InterpolationOrder -> 1]; times = Table[ti, {ti, 0, 20, 0.2}]; interfaceData = Table[hSol @@ {ti, θ}, {ti, times}]; iCs = Table[interfaceCenters[θ, i], {i, interfaceData}]; Qs = Q /@ iCs;  And the animation: ListAnimate@(animation = Table[ Show[ ListPolarPlot[Transpose[{θ, interfaceData[[i]]}], PlotRange -> {{-2, 2}, {-2, 2}}, ImageSize -> Large, PlotStyle -> Black, Joined -> True, Epilog -> {Text["Q = " <> ToString[Qs[[i]]], {0.1, 0.1}, {-1, -1}]} ], RegionPlot[1 <= x^2 + y^2 <= 4, {x, -2, 2}, {y, -2, 2}, PlotStyle -> Opacity[0.1], AspectRatio -> 1 ], ListPolarPlot[{{3 π/2, 0} + {-1, 1} iCs[[i, 1]], {3 π/2, 0} + iCs[[i, 2]]}, PlotStyle -> {{Blue, PointSize -> Large}} ], ListPolarPlot[ Table[{{0, 0}, {{π, 0} + iCs[[i, 1]], {3 π/2, 0} + iCs[[i, 2]]}[[j]]}, {j, Length[iCs[[i]]]} ], PlotStyle -> {{Black, Dashed}}, Joined -> True ] ], {i, Length[times]}]) Export["C:\\Users\\jdbishop\\Desktop\\test.gif", animation]  ## Update I think WhenEvent works basically like it is supposed to, and therefore constitutes a solution that is technically correct, it's just that using it to update$Q$every means that the simulation misses the case when$Q=0$and therefore the dynamics reach equilibrium. I've edited my answer a final time to correct the WhenEvent code to correct the pattern matching failure pointed out by @xzczd. However, you can see the failure to catch$Q=0$in the animation above. So I see @xzczd's answer and think that updating$Q$discretely is never going to be a satisfactory answer when @xzczd has shown how to update it continuously while generating the numerical solution. So that's it from me, I think, but I've learned a lot working on this problem with everyone else. • Is there a reason you changed the initial condition h0 to have the points on the outer and inner walls exactly 1 and 2? They need to be 1.001 and 1.999 respectively or else the PDE will be equal to 0 on the outer and inner walls and will not move there, this gives an incorrect solution. When I try to change the initial condition to the previous correct one (after making the appropriate amends in interfaceCentres I get errors about Part{2,4} not existing. – M Horsley Aug 9 '16 at 7:58 • Also, if I take the simulation further to$t = 60$say, the centers of the interface clearly swap "heights", ie. the flat section on the right moves below the flat section on the right. This must mean that$\psi = \phi$at some point, which in turn must mean$Q = 0$. This is missed by your code as$Q \neq 0$, indeed it does not even go negative as the interfaces are swapped. – M Horsley Aug 9 '16 at 9:39 • I'm still massively confused about why changing the initial condition back to my original one is giving errors. I amend h0 to be my original one, and then amend N@SplitBy[ Transpose[{h\[Theta]t, \[Theta]}], #[[1]] != 1.999 && #[[1]] != 1.001 &][[{2, 4}]] This amend returns the same as your code, which it should do as it throws away the inner/outer pieces and just keeps the flat bits. So given these amends that return the same code as yours, I have no idea why the simulation section is throwing up Part{2,4}` does not exist errors. – M Horsley Aug 9 '16 at 10:36 • @ojlm Sorry! I wasn't sure that the pattern matching approach was working that great, so I used a different approach that clearly turned out to be a rather fragile one. I've edited my answer to incorporate your center-finding function instead, which actually seems to be way faster anyway! The code ins't as general as it was before, but I think it's implemented correctly, now (at least I hope so). The remaining confusion on my side of things is why does the center of the right interface jump so quickly at$t = 0.2 s\$? It doesn't seem like it's finding the correct center position... – Josh Bishop Aug 9 '16 at 17:55
• – Michael E2 Aug 10 '16 at 20:16