如何在 Mathematica 的 NDSolve 中“嵌入”分段答案

【问题标题】：How to "embed" Piecewise in NDSolve in Mathematica如何在 Mathematica 的 NDSolve 中“嵌入”分段
【发布时间】：2011-09-21 18:14:52
【问题描述】：

我正在使用NDSolve 解决非线性偏微分方程。
我希望其中一个变量 (Kvar) 成为函数当前正在解决的时间步长，因此和使用 Piecewise。
Mathematica 生成一条错误消息：

SetDelayed::write: 0.05[t_] 中的实数标记受保护。 >> NDSolve::deqn: 期望的方程或方程列表而不是 $Failed 在第一个参数中......

ReplaceAll::reps: ....

为了便于阅读，我没有包含完整的错误消息。

我的代码如下：

Needs["VectorAnalysis`"]
Needs["DifferentialEquations`InterpolatingFunctionAnatomy`"];
Clear[Eq4, EvapThickFilm, h, S, G, E1, K1, D1, VR, M, R]
Eq4[h_, {S_, G_, E1_, K1_, D1_, VR_, M_, R_}] := \!\(
\*SubscriptBox[\(\[PartialD]\), \(t\)]h\) + 
    Div[-h^3 G Grad[h] + 
      h^3 S Grad[Laplacian[h]] + (VR E1^2 h^3)/(D1 (h + K1)^3)
        Grad[h] + M (h/(1 + h))^2 Grad[h]] + E1/(
    h + K1) + (R/6) D[D[(h^2/(1 + h)), x] h^3, x] == 0;
SetCoordinates[Cartesian[x, y, z]];
EvapThickFilm[S_, G_, E1_, K1_, D1_, VR_, M_, R_] := 
  Eq4[h[x, y, t], {S, G, E1, K1, D1, VR, M, R}];
TraditionalForm[EvapThickFilm[S, G, E1, K1, D1, VR, M, R]];

我尝试在NDSolve 中实现Piecewise 的第二个单元格：

L = 318; TMax = 7.0;
Off[NDSolve::mxsst];
(*Ktemp = Array[0.001+0.001#^2&,13]*)
hSol = h /. NDSolve[{
     (*S,G,E,K,D,VR,M*)

 Kvar[t_] :=  Piecewise[{{0.01, t <= 4}, {0.05, t > 4}}],
 EvapThickFilm[1, 3, 0.1, Kvar[t], 0.01, 0.1, 0, 160],
 h[0, y, t] == h[L, y, t],
 h[x, 0, t] == h[x, L, t],
 (*h[x,y,0] == 1.1+Cos[x] Sin[2y] *)

 h[x, y, 0] == 
  1 + (-0.25 Cos[2 \[Pi] x/L] - 0.25 Sin[2 \[Pi] x/L]) Cos[
     2 \[Pi] y/L]
 },
h,
{x, 0, L},
{y, 0, L},
{t, 0, TMax}
][[1]]

hGrid = InterpolatingFunctionGrid[hSol];

PS：很抱歉，这里的第一个单元格块显示得不太好。而且由于没有足够的“声誉”，我无法发布图片。

使用NDSolve 单元格块时出现错误消息。

【问题讨论】：

从 NDSolve 参数列表中取出 Kvar[t_] 的定义。然后它会给出一个解决方案。可能需要调整 NDSolve 设置以获得更好的错误估计并成功运行到 TMax（似乎在 2.36 左右退出）。
谢谢！这当然有帮助。就像你建议的那样，接下来我将尝试使用 NDSolve 设置。有什么建议吗？
最好不要使用以大写字母开头的变量名称，因为这可能会导致与内置名称（均以大写字母开头）发生冲突。 TraditionalForm[EvapThickFilm[S, G, E1, K1, D1, VR, M, R]]; 这行没有意义。没有分配和 ;抑制任何输出。因此没有净效应。
我取消“TraditionalForm[...]”的唯一原因是因为输出太长而且不好看。

标签： wolfram-mathematica piecewise

【解决方案1】：

在NDSolve 中的一组方程之外定义函数Kvar，如

Off[NDSolve::mxsst];
(*Ktemp=Array[0.001+0.001#^2&,13]*)
Kvar[t_] := Piecewise[{{0.01, t <= 4}, {0.05, t > 4}}];
hSol = ...

并将其从NDSolve 的列表中删除，使其以NDSolve[{(*S,G,E,K,D,VR,M*)EvapThickFilm[... 开头，这样就可以工作了。它会给出警告，但这些警告与方程中可能存在的奇点有关。

此外，您的原始错误表明您的Kvar 被分配了0.05 的值。因此，请在第二个单元格中的任何其他内容之前添加 Clear[Kvar]。

【讨论】：