首先,LZ对分段的初始条件的定义方法显然是错误的,正确的做法恰恰是,就按你图上的来:
DiffusionEquation = {D[\[Rho][r, t], t] - (d D[D[\[Rho][r, t], r] r, r])/r ==
0, \[Rho][r, 0] == Piecewise[{{c, r < R0}}]};
然后,目前DSolve在处理无限域问题时的语法并不是在自变量定义域的位置写Infinity,而是,就直接空着。这部分内容参看DSolve的自带帮助的“范围”部分,这里不展开讲。因为——上面这两点问题相较于这个问题的真正麻烦之处根本就是微不足道。
这个问题的真正麻烦之处在于:
1. DSolve在经历了10.3版的增强之后,似乎依旧不能直接求解这个问题。
2. 极坐标无限域下的Laplacian其实可以使用汉克尔变换消掉微分,但是遗憾的是目前Mathematica的HankelTransform性能还很糟,当然了,这个还不算太难变通。
3. 最终我们所得的变换后的解,不知道是太复杂了怎么的,总之,倒不回来。当然了,我们依旧得到了一个含了积分号的解析解。
好的,我们开始求解。首先利用拉普拉斯变换消掉关于t的导数项:
DiffusionEquation = {D[\[Rho][r, t], t] - (d D[D[\[Rho][r, t], r] r, r])/r ==
0, \[Rho][r, 0] == Piecewise[{{c, r < R0}}]};
teqn = LaplaceTransform[{DiffusionEquation[[1]]}, t, s] /.
Rule @@ DiffusionEquation[[2]] /.
HoldPattern@LaplaceTransform[a_, __] :> a /. \[Rho] -> (rho[#] &)
(*
这里做了一些字母替换,主要是为了方便阅读。大家记住rho[r]代表\[Rho][r, t]对t的拉普拉斯变换就行了。然后呢,假如HankelTransform具备和LaplaceTransform相当的性能的话,我们只要使用类似上面的步骤,就可以消掉关于r的导数了,但是,遗憾的是:
*)
-d HankelTransform[rho'[r]/r + rho''[r], r, w]
HankelTransform[-(rho'[r]/r + rho''[r]), r, w]
(*嗯你没看错,只是加了个负号,HankelTransform就跪了。更囧的是:*)
Assuming[{R0 > 0}, HankelTransform[Piecewise[{{c, r < R0}}, 0] , r, w]]
(*返回未计算表达式*)
(*然而这个算式和自带帮助里的例子其实是差不多的:*)
(*并且Integrate其实也是会算这个积分的:*)
term = Assuming[{R0 > 0},
Integrate[r Piecewise[{{c, r < R0}}, 0] BesselJ[0, r w], {r, 0, Infinity}]]
(* (c R0 BesselJ[1, R0 w])/w *)
(*
总之呢,经历以上步骤我们总算是手动拼出了Hankel变换后的方程。我们把它解出来:
*)
ttsol = p /. First@Solve[-term + s p + d w^2 p == 0, p]
(* (c R0 BesselJ[1, R0 w])/(w (s + d w^2)) *)
(*然后拉普拉斯逆变换:*)
tsol = InverseLaplaceTransform[ttsol, s, t]
(*
(c E^(-d t w^2) R0 BesselJ[1, R0 w])/w
*)
(*
现在,我们来到了最后的难关,汉克尔逆变换。很遗憾,这个符号解我们算不出来:
*)
InverseHankelTransform[tsol, w, r]
Block[{d = 1, c = 1, R0 = 1, t = 1}, InverseHankelTransform[tsol, w, r]]
Block[{d = 1, c = 1, R0 = 1, t = 1}, Integrate[w tsol BesselJ[0, w r], {w, 0, Infinity}]]
(*以上方法全部失效*)
(*不过我们总算是得到了一个含有积分号的解析解:*)
(*这个解析解当然也是可以用的。我们和数值解对比一下:*)
asol[t_, r_, c_, d_, R0_?NumericQ] :=
NIntegrate[w (c E^(-d t w^2) R0 BesselJ[1, R0 w])/w BesselJ[0, w r], {w, 0, Infinity}]
aplot = Plot[asol[1, r, 1, 1, 1], {r, 0, 3}, PlotRange -> All]; // AbsoluteTiming
(* {6.04077, Null} *)
nsol = Block[{d = 1, c = 1, R0 = 1},
NDSolveValue[DiffusionEquation, \[Rho], {t, 0, 2}, {r, 0, 10},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.01}}]]
(*注意,当有限元方法被用于NDSolve的空间离散时,如果不明确指定边界,NDSolve会自动使用等于0的NeumannValue做为边界——我们在此需要的就是这种边界。*)
(* 对比一下: *)
aplot~Show~Plot[nsol[r, 1], {r, -3, 3}, PlotStyle -> {Red, Dashed}]
DiffusionEquation = {D[\[Rho][r, t], t] - (d D[D[\[Rho][r, t], r] r, r])/r ==
0, \[Rho][r, 0] == Piecewise[{{c, r < R0}}]};
然后,目前DSolve在处理无限域问题时的语法并不是在自变量定义域的位置写Infinity,而是,就直接空着。这部分内容参看DSolve的自带帮助的“范围”部分,这里不展开讲。因为——上面这两点问题相较于这个问题的真正麻烦之处根本就是微不足道。
这个问题的真正麻烦之处在于:
1. DSolve在经历了10.3版的增强之后,似乎依旧不能直接求解这个问题。
2. 极坐标无限域下的Laplacian其实可以使用汉克尔变换消掉微分,但是遗憾的是目前Mathematica的HankelTransform性能还很糟,当然了,这个还不算太难变通。
3. 最终我们所得的变换后的解,不知道是太复杂了怎么的,总之,倒不回来。当然了,我们依旧得到了一个含了积分号的解析解。
好的,我们开始求解。首先利用拉普拉斯变换消掉关于t的导数项:
DiffusionEquation = {D[\[Rho][r, t], t] - (d D[D[\[Rho][r, t], r] r, r])/r ==
0, \[Rho][r, 0] == Piecewise[{{c, r < R0}}]};
teqn = LaplaceTransform[{DiffusionEquation[[1]]}, t, s] /.
Rule @@ DiffusionEquation[[2]] /.
HoldPattern@LaplaceTransform[a_, __] :> a /. \[Rho] -> (rho[#] &)
(*
这里做了一些字母替换,主要是为了方便阅读。大家记住rho[r]代表\[Rho][r, t]对t的拉普拉斯变换就行了。然后呢,假如HankelTransform具备和LaplaceTransform相当的性能的话,我们只要使用类似上面的步骤,就可以消掉关于r的导数了,但是,遗憾的是:
*)
-d HankelTransform[rho'[r]/r + rho''[r], r, w]
HankelTransform[-(rho'[r]/r + rho''[r]), r, w]
(*嗯你没看错,只是加了个负号,HankelTransform就跪了。更囧的是:*)
Assuming[{R0 > 0}, HankelTransform[Piecewise[{{c, r < R0}}, 0] , r, w]]
(*返回未计算表达式*)
(*然而这个算式和自带帮助里的例子其实是差不多的:*)
(*并且Integrate其实也是会算这个积分的:*)
term = Assuming[{R0 > 0},
Integrate[r Piecewise[{{c, r < R0}}, 0] BesselJ[0, r w], {r, 0, Infinity}]]
(* (c R0 BesselJ[1, R0 w])/w *)
(*
总之呢,经历以上步骤我们总算是手动拼出了Hankel变换后的方程。我们把它解出来:
*)
ttsol = p /. First@Solve[-term + s p + d w^2 p == 0, p]
(* (c R0 BesselJ[1, R0 w])/(w (s + d w^2)) *)
(*然后拉普拉斯逆变换:*)
tsol = InverseLaplaceTransform[ttsol, s, t]
(*
(c E^(-d t w^2) R0 BesselJ[1, R0 w])/w
*)
(*
现在,我们来到了最后的难关,汉克尔逆变换。很遗憾,这个符号解我们算不出来:
*)
InverseHankelTransform[tsol, w, r]
Block[{d = 1, c = 1, R0 = 1, t = 1}, InverseHankelTransform[tsol, w, r]]
Block[{d = 1, c = 1, R0 = 1, t = 1}, Integrate[w tsol BesselJ[0, w r], {w, 0, Infinity}]]
(*以上方法全部失效*)
(*不过我们总算是得到了一个含有积分号的解析解:*)
(*这个解析解当然也是可以用的。我们和数值解对比一下:*)
asol[t_, r_, c_, d_, R0_?NumericQ] :=
NIntegrate[w (c E^(-d t w^2) R0 BesselJ[1, R0 w])/w BesselJ[0, w r], {w, 0, Infinity}]
aplot = Plot[asol[1, r, 1, 1, 1], {r, 0, 3}, PlotRange -> All]; // AbsoluteTiming
(* {6.04077, Null} *)
nsol = Block[{d = 1, c = 1, R0 = 1},
NDSolveValue[DiffusionEquation, \[Rho], {t, 0, 2}, {r, 0, 10},
Method -> {"MethodOfLines",
"SpatialDiscretization" -> {FiniteElement, MeshOptions -> MaxCellMeasure -> 0.01}}]]
(*注意,当有限元方法被用于NDSolve的空间离散时,如果不明确指定边界,NDSolve会自动使用等于0的NeumannValue做为边界——我们在此需要的就是这种边界。*)
(* 对比一下: *)
aplot~Show~Plot[nsol[r, 1], {r, -3, 3}, PlotStyle -> {Red, Dashed}]