Question

用mathematica求解具有边界条件的连续性方程（偏微分方程）,要可执行的代码

Answer 1

……上面那个答案明明什么都没说，为什么你采纳了啊？

算了不管了，总之这里是我的答案。这个问题的关键在于DiracDelta的处理。DiracDelta可以做近似处理，比如说使用：

eq = With[{n = n[t, x]}, D[n, t] == D[n, x, x]];

bc = {n[t, 0] == 0, D[n[t, x], x] == 0 /. x -> 1};

dirac[r_, a_] := Sqrt[a/Pi] Exp[-a r^2]

ic = n[0, x] == 2 n0 dirac[x - 1, 100];(*注意这个 2 是必要的，因为你的DiracDelta设在了端点上*)

nsol = NDSolveValue[{eq, ic, bc} /. n0 -> 1, n, {t, 0, 1/10}, {x, 0, 1}];

np = Plot3D[nsol[t, x], {t, 0, 1/10}, {x, 0, 1}, PlotRange -> All, PlotStyle -> Red];

效果就已经不错了。

不过，对于你这个问题，还存在更准确的处理方法，因为它的解析解是可求的，只要用上LaplaceTransform:

    eq = With[{n = n[t, x]}, D[n, t] == D[n, x, x]];

    bc = {n[t, 0] == 0, D[n[t, x], x] == 0 /. x -> 1};

    ic = n[0, x] == n0 DiracDelta[x - 1];

    

    teqn = 

     LaplaceTransform[{eq, bc}, t, s] /. Rule @@ ic /. 

       HoldPattern@LaplaceTransform[a_, __] :> a /. n -> (n@#2 &)

    (* {-n0 DiracDelta[-1 + x] + s n[x] == (n^′′)[x], {n[0] == 0, 

      Derivative[1][n][1] == 0}} *)

    

    tsol = n[x] /. First@DSolve[teqn, n[x], x]

    (* (1/(2 (1 + E^(2 Sqrt[s])) Sqrt[s]))E^(-Sqrt[s] - 

      Sqrt[s] x) n0 (-2 E^(2 Sqrt[s]) + 2 E^(2 Sqrt[s] + 2 Sqrt[s] x) + 

       E^(2 Sqrt[s]) HeavisideTheta[-1 + x] + E^(4 Sqrt[s]) HeavisideTheta[-1 + x] - 

       E^(2 Sqrt[s] x) HeavisideTheta[-1 + x] - 

       E^(2 Sqrt[s] + 2 Sqrt[s] x) HeavisideTheta[-1 + x]) *)

tsol 是经过了拉普拉斯变换的这个问题的解析解，最后一步就是倒回去了，不过呢它的符号求解似乎不可能，于是我们如果要数值解的话，就需要数值拉普拉斯反变换。这里我将使用这个程序包：library.wolfram.com/infocenter/MathSource/5026

这里以n0等于1为例：

tsolfunc = Function[{s, x}, #] &[tsol /. n0 -> 1]

(*注意FT要从上面的程序包里调出来。*)

solgenerator[t_] := 

 solgenerator[t] = Module[{x}, Compile[#, #2] & @@ {x, FT[tsolfunc[#, x] &, t]}]

Plot3D[solgenerator[t][x], {t, 0, 1/10}, {x, 0, 1}, PlotRange -> All]

最后，这个帖子里是个类似的问题。其中有一些针对DiracDelta的更详细说明，有兴趣可以读读：http://tieba.baidu.com/p/4126799131

Answer 2

Please take a look my example. Besides, check your help document before posting your question.
Using D to take derivatives
pde = D[y[x, t], t] + 2 D[y[x, t], x] == 0
Use DSolve to solve the equation and store the solution as soln.
soln = DSolve[pde, y[x, t], {x, t}]
To use the solution as a function, say f[x,t], use /. (the short form of ReplaceAll) and [[...]]
f[x_, t_] = y[x, t] /. soln[[1]]