Skip to main content

differential equations - How to use NDSolve with moving boundary conditions?


So I am trying to solve the movement in space and time of a spreading gravity current. The interface satisfies the following PDE:



$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) $.


The nose, $h(x,t) = 0$, of the current can move forward in space. If we say the initial condition is $h(x,0) = 1-x$, then the nose, $x_N$ is initially at $x = 1$. The base of the current can be fixed at $h(0,t) = 1$. So I have one IC and one BC, I need another BC to close the system. Now, as time moves forward the nose (fixed at h = 0) propagates according to the expression:


$\int_0^{x_N} h \, dx = \frac{1}{2} + t$


or kinematically $\dot{x_N} = -h^2 h_x$.


Clearly the second boundary condition in $x$ needs to come from the nose condition, but I have no idea how to input an integral boundary condition or the kinematic condition. Secondly, I have no idea how to deal with xmax in NDSolve as clearly my xmax is moving forward with each time step.


Here is my attempt, which gives a "There are fewer dependent variables error"


NDSolve[{D[h[x, t], t] == D[h[x, t]^3 D[h[x, t], x], x], h[0, t] == 1,
h[x, 0] == 1 - x (D[h[x, t], t] /. x -> 1) == -h[x, t]^2 (D[h[x, t], x] /.
x -> 1)}, h, {x, 0, 1}, {t, 0, 1}]


EDIT


So I solved this problem by amending the code in this Stefan problem post. I modified the statement of the problem slightly to work more closely to @ybeltukov excellent solution, and actually solved the following system:


$ \frac{\partial h}{\partial t} = \frac{\partial}{\partial x}\left(h^3 \frac{\partial h}{\partial x}\right) \\ \dot{s} = -h^2 h_x \\ s(0) = 0 \\ h(x,0) = \begin{cases}1 \quad \text{for } x = 0 \\ 0 \quad \text{otherwise}\end{cases} \\ h(s(t),t) = 0 \\ h(0,t) = 1 $.


I made the same change to a normalised variable for my PDE and then amended finite difference code:


n = 100;
\[Delta]\[Xi] = 1./n;

ClearAll[dv, t];
dv[v_List] :=
With[{s = First@v, u = Rest@v},

With[{ds = u[[-1]]^3/(s \[Delta]\[Xi]), \[Xi] = N@Range[n - 1]/n,
d1 = ListCorrelate[{-0.5, 0, 0.5}/\[Delta]\[Xi], #] &,
d2 = ListCorrelate[{1, -2, 1}/\[Delta]\[Xi]^2, #] &},
Prepend[3 u^2 d1[#]^2/s^2 + u^3 d2[#]/s^2 + \[Xi] ds d1[#]/s &@
Join[{1}, u, {0.}], ds]]];
s0 = 0.001;
v0 = Flatten@Prepend[ConstantArray[0.001, n - 2], {s0, 1.}];
sol = NDSolve[{v'[t] == dv[v[t]], v[0] == v0}, v, {t, 0, 50}][[1, 1,
2]];


Returning from the normalised variable is identical code.


You can then make a nice time series with this code:


ListAnimate[
Table[Plot[u[t, x], {x, 0, 10}, PlotRange -> {{0, 10}, {0, 3}}], {t,
0, 50, 0.1}]]

Which returns a spreading current as desired.


Spreading viscous gravity current




Comments

Popular posts from this blog

front end - keyboard shortcut to invoke Insert new matrix

I frequently need to type in some matrices, and the menu command Insert > Table/Matrix > New... allows matrices with lines drawn between columns and rows, which is very helpful. I would like to make a keyboard shortcut for it, but cannot find the relevant frontend token command (4209405) for it. Since the FullForm[] and InputForm[] of matrices with lines drawn between rows and columns is the same as those without lines, it's hard to do this via 3rd party system-wide text expanders (e.g. autohotkey or atext on mac). How does one assign a keyboard shortcut for the menu item Insert > Table/Matrix > New... , preferably using only mathematica? Thanks! Answer In the MenuSetup.tr (for linux located in the $InstallationDirectory/SystemFiles/FrontEnd/TextResources/X/ directory), I changed the line MenuItem["&New...", "CreateGridBoxDialog"] to read MenuItem["&New...", "CreateGridBoxDialog", MenuKey["m", Modifiers-...

How to thread a list

I have data in format data = {{a1, a2}, {b1, b2}, {c1, c2}, {d1, d2}} Tableform: I want to thread it to : tdata = {{{a1, b1}, {a2, b2}}, {{a1, c1}, {a2, c2}}, {{a1, d1}, {a2, d2}}} Tableform: And I would like to do better then pseudofunction[n_] := Transpose[{data2[[1]], data2[[n]]}]; SetAttributes[pseudofunction, Listable]; Range[2, 4] // pseudofunction Here is my benchmark data, where data3 is normal sample of real data. data3 = Drop[ExcelWorkBook[[Column1 ;; Column4]], None, 1]; data2 = {a #, b #, c #, d #} & /@ Range[1, 10^5]; data = RandomReal[{0, 1}, {10^6, 4}]; Here is my benchmark code kptnw[list_] := Transpose[{Table[First@#, {Length@# - 1}], Rest@#}, {3, 1, 2}] &@list kptnw2[list_] := Transpose[{ConstantArray[First@#, Length@# - 1], Rest@#}, {3, 1, 2}] &@list OleksandrR[list_] := Flatten[Outer[List, List@First[list], Rest[list], 1], {{2}, {1, 4}}] paradox2[list_] := Partition[Riffle[list[[1]], #], 2] & /@ Drop[list, 1] RM[list_] := FoldList[Transpose[{First@li...

plotting - Magnifying Glass on a Plot

Although there is a trick in TEX magnifying glass but I want to know is there any function to magnifying glass on a plot with Mathematica ? For example for a function as Sin[x] and at x=Pi/6 Below, this is just a picture desired from the cited site. the image got huge unfortunately I don't know how can I change the size of an image here! Answer Insetting a magnified part of the original Plot A) by adding a new Plot of the specified range xPos = Pi/6; range = 0.2; f = Sin; xyMinMax = {{xPos - range, xPos + range}, {f[xPos] - range*GoldenRatio^-1, f[xPos] + range*GoldenRatio^-1}}; Plot[f[x], {x, 0, 5}, Epilog -> {Transparent, EdgeForm[Thick], Rectangle[Sequence @@ Transpose[xyMinMax]], Inset[Plot[f[x], {x, xPos - range, xPos + range}, Frame -> True, Axes -> False, PlotRange -> xyMinMax, ImageSize -> 270], {4., 0.5}]}, ImageSize -> 700] B) by adding a new Plot within a Circle mf = RegionMember[Disk[{xPos, f[xPos]}, {range, range/GoldenRatio}]] Show...