Skip to main content

numerical integration - Heat equation with nonlinear boundary condition involving time-derivative


The governing equation is shown as follows:


enter image description here


I first try to employ the NDSolve, but it seems that Mathematica can not handle the fourth boundary condition.


Therefore, I rewrite the code with finite difference method, but still have problem


Mr = 100; Mz = 10; Lr = 10; Lz = 1;
γ = 0.01;
α = 1;

κ = 1;
dr = Lr/Mr;
dz = Lz/dz;
dt = 0.1;
r[i_] := i*dr;
z[j_] := j*dz;


(*i.c.*)
u[i_, j_, 0] := 0;

(*b.c.*)
u[1, j_, k_] := u[2, j, k] + 2;
u[Mr, j_, k_] := u[Mr - 1, j, k];
u[i_, 1, k_] := u[i, 2, k];
u[i_, Mz, k_] :=
If[k == 1, 0,
1/dz^2 (dt dz γ u[i, -1 + Mz, -1 + k] +
dt α γ u[i, -1 + Mz, -1 + k]^2 +
dz^2 u[i, Mz, -1 + k] - dt dz γ u[i, Mz, -1 + k] -
2 dt α γ u[i, -1 + Mz, -1 + k] u[i, Mz, -1 + k] +

dt α γ u[i, Mz, -1 + k]^2)];

u[i_, j_, k_] :=
u[i, j, k] =
1/(dr^2 dz^2 i) (dt dz^2 i u[-1 + i, j, -1 + k] +
dr^2 dt i u[i, -1 + j, -1 + k] - dt dz^2 u[i, j, -1 + k] -
2 dr^2 dt i u[i, j, -1 + k] + dr^2 dz^2 i u[i, j, -1 + k] -
2 dt dz^2 i u[i, j, -1 + k] + dr^2 dt i u[i, 1 + j, -1 + k] +
dt dz^2 u[1 + i, j, -1 + k] + dt dz^2 i u[1 + i, j, -1 + k]);


u[1, Mz, 1]

However, I receive the following error message



$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>


$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>


$RecursionLimit::reclim: Recursion depth of 1024 exceeded. >>


General::stop: Further output of $RecursionLimit::reclim will be suppressed during this calculation. >>



The problem may be induced by the imposion of nonlinear and time-derivative and nonlinear boundary condition. Any idea about how to solve my problem? This boundary condition has bothered me for long time.





I re code the provided code as follows


Clear[fdd, pdetoode, tooderule]
fdd[{}, grid_, value_, order_] := value;
fdd[a__] := NDSolve`FiniteDifferenceDerivative@a;

pdetoode[funcvalue_List, rest__] :=
pdetoode[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]],
rest];
pdetoode[{func__}[var__], rest__] :=

pdetoode[Alternatives[func][var], rest];
pdetoode[rest__, grid_?VectorQ, o_Integer] :=
pdetoode[rest, {grid}, o];

pdetoode[func_[var__], time_, {grid : {__} ..}, o_Integer] :=
With[{pos = Position[{var}, time][[1, 1]]},
With[{bound = #[[{1, -1}]] & /@ {grid},
pat = Repeated[_, {pos - 1}],
spacevar = Alternatives @@ Delete[{var}, pos]},
With[{coordtoindex =

Function[coord,
MapThread[
Piecewise[{{1, # === #2[[1]]}, {-1, # === #2[[-1]]}},
All] &, {coord, bound}]]},
tooderule@
Flatten@{((u : func) |
Derivative[dx1 : pat, dt_, dx2___][(u : func)])[x1 : pat,
t_, x2___] :> (Sow@coordtoindex@{x1, x2};

fdd[{dx1, dx2}, {grid},

Outer[Derivative[dt][u@##]@t &, grid],
"DifferenceOrder" -> o]),
inde : spacevar :>
With[{i = Position[spacevar, inde][[1, 1]]},
Outer[Slot@i &, grid]]}]]];

tooderule[rule_][pde_List] := tooderule[rule] /@ pde;
tooderule[rule_]@Equal[a_, b_] :=
Equal[tooderule[rule][a - b], 0] //.
eqn : HoldPattern@Equal[_, _] :> Thread@eqn;

tooderule[rule_][expr_] := #[[Sequence @@ #2[[1, 1]]]] & @@
Reap[expr /. rule]

Clear@pdetoae;
pdetoae[funcvalue_List, rest__] :=
pdetoae[(Alternatives @@ Head /@ funcvalue) @@ funcvalue[[1]], rest];
pdetoae[{func__}[var__], rest__] :=
pdetoae[Alternatives[func][var], rest];

pdetoae[func_[var__], rest__] :=

Module[{t},
Function[
pde, #[pde /. {Derivative[d__][u : func][inde__] :>
Derivative[d, 0][u][inde, t], (u : func)[inde__] :>
u[inde, t]}] /. (u : func)[i__][t] :> u[i]] &@
pdetoode[func[var, t], t, rest]];
\[Gamma] = 100;
\[Alpha] = 0;
\[Kappa] = 1;
R = 10;

Z = 1;
eps = 10^-1;
tend = 1000;

eq = With[{u = u[r, Sqrt[\[Kappa]] z, t]},
Laplacian[u, {r, th, z}, "Cylindrical"] == D[u, t] /.
Sqrt[\[Kappa]] z -> z];

ic = u[r, z, 0] == 0;
bc = With[{u = u[r, z, t]}, {D[u, r] == -2/r /. r -> eps,

u == 0 /. r -> R, D[u, z] == 0 /. z -> 0, D[u, z] == 0 /. z -> Z}]

domain@r = {eps, R};
domain@z = {0, Z};
points@r = 50;
points@z = 50;
difforder = 2;
(grid@# = Array[# &, points@#, domain@#]) & /@ {r, z};

(*Definition of pdetoode isn't included in this post,please find it \

in the link above.*)
ptoofunc = pdetoode[u[r, z, t], t, grid /@ {r, z}, difforder];

delbothside = #[[2 ;; -2]] &;

ode = delbothside /@ delbothside@ptoofunc@eq;

odeic = delbothside /@ delbothside@ptoofunc@ic;
odebc = MapAt[delbothside, ptoofunc@bc, {{1}, {2}}];


sollst = NDSolveValue[{ode, odeic, odebc},
Outer[u, grid@r, grid@z] // Flatten, {t, 0, tend}];

sol = ListInterpolation[
Partition[Developer`ToPackedArray@#["ValuesOnGrid"] & /@ #,
points@z], {grid@r, grid@z, #[[1]]["Coordinates"][[1]]}] &@
sollst; // AbsoluteTiming

Manipulate[
Plot3D[sol[r, z, t], {r, eps, R}, {z, 0, Z},

PlotRange -> {-1, 10}], {t, 0, tend}]

The equation indicate that the upper and under boundary conditions are subjected to no-flow condition. And it can easily be solved and obtain an analytical solution.


To compare the analytical solution (or semi-analytical), the solution can be expressed as


Vi[n_, i_] := 
Vi[n, i] = (-1)^(i + n/2) Sum[
k^(n/2) (2 k)! /( (n/2 - k)! k! (k - 1)! (i - k)! (2 k -
i)! ), { k, Floor[ (i + 1)/2 ], Min[ i, n/2] } ] // N;
Stehfest[F_, s_, t_, n_: 16] :=
If[n > 16, Message[Stehfest::optimalterms, n];

If[ OddQ[n], Message[Stehfest::odd, n];
"Enter an even number of terms",
If[n > 32, Message[Stehfest::terms, n];
" Try a smaller value for n. Maximum allowable n is 32 ",
Log[2]/t Sum[ Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]],
If[ OddQ[n], Message[Stehfest::odd, n];
"Enter an even number of terms",
If[n > 32, Message[Stehfest::terms, n];
" Try a smaller value for n. Maximum allowable n is 32.",
Log[2]/t Sum[

Vi[n, i]*F /. s -> i Log[2]/t , {i, 1, n} ] ]]] // N;
s0[r_, z_, t_] :=
NIntegrate[
Re[Stehfest[(
2 E^(-((Sqrt[a^2 + p] z)/
Sqrt[\[Kappa]])) (-E^((Sqrt[a^2 + p]/Sqrt[\[Kappa]]))
p \[Gamma] -
E^((Sqrt[a^2 + p] (1 + 2 z))/Sqrt[\[Kappa]]) p \[Gamma] +
E^((Sqrt[a^2 + p] z)/
Sqrt[\[Kappa]]) (p \[Gamma] - Sqrt[a^2 + p] Sqrt[\[Kappa]]) +

E^((Sqrt[a^2 + p] (2 + z))/
Sqrt[\[Kappa]]) (p \[Gamma] +
Sqrt[a^2 + p] Sqrt[\[Kappa]])))/(
p (a^2 +
p) ((1 + E^((2 Sqrt[a^2 + p])/
Sqrt[\[Kappa]])) p \[Gamma] + (-1 + E^((2 Sqrt[a^2 + p])/
Sqrt[\[Kappa]])) Sqrt[a^2 + p] Sqrt[\[Kappa]])), p, t, 6]]*
BesselJ[0, a (r)]*a, {a, 0, \[Infinity]},
Method -> {"LevinRule", "LevinFunctions" -> {BesselJ}},
MaxRecursion -> 40];


with \[Gamma] = 0; I however compare the numerical results and analytical one. It seems that the the numerical solution obtain a wrong result. The code and the plot of comparison are shown as follows


\[Gamma] = 0; TT1 = Table[{t, s0[0.9, 0.5, t]}, {t, 0.1, 1, 0.1}];
TT2 = Table[{t, s0[0.9, 0.5, t]}, {t, 1, 10, 1}];
TT3 = Table[{t, s0[0.9, 0.5, t]}, {t, 10, 100, 10}];
TT4 = Table[{t, s0[0.9, 0.5, t]}, {t, 100, 1000, 100}];
C0 = Join[TT1, TT2, TT3, TT4];
TT1 = Table[{t, sol[1, 0.5, t]}, {t, 0.1, 1, 0.1}];
TT2 = Table[{t, sol[1, 0.5, t]}, {t, 1, 10, 1}];
TT3 = Table[{t, sol[1, 0.5, t]}, {t, 10, 100, 10}];

TT4 = Table[{t, sol[1, 0.5, t]}, {t, 100, 1000, 100}];
Cn = Join[TT1, TT2, TT3, TT4];

enter image description here


The figure is u versus t and blue and yellow lines represent the value predicted by numerical and analytical solutions, respectively. I an trying to figure out what happen to the numerical method.




Comments

Popular posts from this blog

mathematical optimization - Minimizing using indices, error: Part::pkspec1: The expression cannot be used as a part specification

I want to use Minimize where the variables to minimize are indices pointing into an array. Here a MWE that hopefully shows what my problem is. vars = u@# & /@ Range[3]; cons = Flatten@ { Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; Minimize[{Total@((vec1[[#]] - vec2[[u[#]]])^2 & /@ Range[1, 3]), cons}, vars, Integers] The error I get: Part::pkspec1: The expression u[1] cannot be used as a part specification. >> Answer Ok, it seems that one can get around Mathematica trying to evaluate vec2[[u[1]]] too early by using the function Indexed[vec2,u[1]] . The working MWE would then look like the following: vars = u@# & /@ Range[3]; cons = Flatten@{ Table[(u[j] != #) & /@ vars[[j + 1 ;; -1]], {j, 1, 3 - 1}], 1 vec1 = {1, 2, 3}; vec2 = {1, 2, 3}; NMinimize[ {Total@((vec1[[#]] - Indexed[vec2, u[#]])^2 & /@ R...

functions - Get leading series expansion term?

Given a function f[x] , I would like to have a function leadingSeries that returns just the leading term in the series around x=0 . For example: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x)] x and leadingSeries[(1/x + 2 + (1 - 1/x^3)/4)/(4 + x)] -(1/(16 x^3)) Is there such a function in Mathematica? Or maybe one can implement it efficiently? EDIT I finally went with the following implementation, based on Carl Woll 's answer: lds[ex_,x_]:=( (ex/.x->(x+O[x]^2))/.SeriesData[U_,Z_,L_List,Mi_,Ma_,De_]:>SeriesData[U,Z,{L[[1]]},Mi,Mi+1,De]//Quiet//Normal) The advantage is, that this one also properly works with functions whose leading term is a constant: lds[Exp[x],x] 1 Answer Update 1 Updated to eliminate SeriesData and to not return additional terms Perhaps you could use: leadingSeries[expr_, x_] := Normal[expr /. x->(x+O[x]^2) /. a_List :> Take[a, 1]] Then for your examples: leadingSeries[(1/x + 2)/(4 + 1/x^2 + x), x] leadingSeries[Exp[x], x] leadingSeries[(1/x + 2 + (1 - 1/x...

What is and isn't a valid variable specification for Manipulate?

I have an expression whose terms have arguments (representing subscripts), like this: myExpr = A[0] + V[1,T] I would like to put it inside a Manipulate to see its value as I move around the parameters. (The goal is eventually to plot it wrt one of the variables inside.) However, Mathematica complains when I set V[1,T] as a manipulated variable: Manipulate[Evaluate[myExpr], {A[0], 0, 1}, {V[1, T], 0, 1}] (*Manipulate::vsform: Manipulate argument {V[1,T],0,1} does not have the correct form for a variable specification. >> *) As a workaround, if I get rid of the symbol T inside the argument, it works fine: Manipulate[ Evaluate[myExpr /. T -> 15], {A[0], 0, 1}, {V[1, 15], 0, 1}] Why this behavior? Can anyone point me to the documentation that says what counts as a valid variable? And is there a way to get Manpiulate to accept an expression with a symbolic argument as a variable? Investigations I've done so far: I tried using variableQ from this answer , but it says V[1...