如何绘制此 PDE 的解?

how to plot the solution to this PDE?

Maple 为这个偏微分方程生成了一个奇怪的解形式。我很难制定解决方案。

解是无穷级数。我将项数设置为 20,然后将时间设置为在 t=2 秒处绘制解决方案。然后现在想为 x=0..1 绘制解决方案。但是情节是空洞的。

当我对解决方案进行采样并使用 listplot 时,我得到了正确的解决方案图。

这是MWE

restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x)+x;
bc:=u(0,t)=0,u(1,t)=0;
ic:=u(x,0)=x*(1-x);
sol:=pdsolve({pde,ic,bc},u(x,t)):
sol:=value(sol);

现在将术语数设置为 20 并设置 t=2

sol2:=subs(t=2,sol):
sol2:=subs(infinity=20,sol2);

以上就是我要的剧情了

plot(rhs(sol2),x=0..1);

我得到空图

所以不得不手动采样并使用 listplot

f:=x->rhs(sol2);
data:=[seq([x,f(x)],x=0..1,.01)]:
plots:-listplot(data);

当我将其与 Mathematica 的结果进行比较时,解决方案看起来是正确的。但是 Mathematica 结果更简单,因为它在总和中没有那些积分。

pde=D[u[x,t],t]==D[u[x,t],{x,2}]+x;
bc={u[0,t]==0,u[1,t]==0};
ic=u[x,0]==x(1-x);
DSolve[{pde,ic,bc},u[x,t],x,t];
%/.K[1]->n;
%/.Infinity->20;
%/.t->2;

剧情是

问题是:如何在不手动采样的情况下绘制 Maple 解决方案?

简短的回答似乎是 Maple 2017.3 中的回归。

对我来说,您的代码可以直接在 Maple 2017.2 和 Maple 2016.2 中运行(没有任何未计算的积分)。我将针对回归提交错误报告。

[编辑]让我知道这四种方法中的任何一种是否适用于您的版本(大概是 Maple 2017.3)。

restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x)+x;
bc:=u(0,t)=0,u(1,t)=0;
ic:=u(x,0)=x*(1-x);
sol:=pdsolve({pde,ic,bc},u(x,t)):
sol:=value(sol);

sol5:=value(combine(subs([sum=Sum,t=2,infinity=20],sol))):
plot(rhs(sol5),x=0..1);

sol4:=combine(subs([sum=Sum,t=2,infinity=20],sol)):
(UseHardwareFloats,oldUHF):=false,UseHardwareFloats:
plot(rhs(sol4),x=0..1);
UseHardwareFloats:=oldUHF: # re-instate

sol2:=subs([sum=Sum,int=Int,t=2],sol):
# Switch integration and summation in second summand of rhs(sol).
sol3:=subsop(2=Sum(int(op([2,1,1],rhs(sol2)),op([2,2],rhs(sol2))),
                   op([2,1,2],rhs(sol2))),rhs(sol2)):
# Rename dummy index and combine summations.
sol3:=Sum(subs(n1=n,op([1,1],sol3))+op([2,1],sol3),
          subs(n1=n,op([1,2],sol3))):
# Curtail to first 20 terms.
sol3:=lhs(sol2)=subs(infinity=20,simplify(sol3));
plot(rhs(sol3),x=0..1);

F:=unapply(subs([Sum='add'],rhs(sol3)),x):
plot(F,0..1);

[edited] 这是另一种方式,在 64 位的 Maple 2017.3 中为我工作Linux。

它可以快速生成图,并且不涉及在 20 项时削减任何金额。请注意,它不会执行 sol:=value(sol); 之前的步骤,因为它会在使用 value 命中任何 Sum 之前激活 int 而不是 Int。它还使用与绘图范围相对应的 x 假设。

restart;
pde:=diff(u(x,t),t)=diff(u(x,t),x)+x:
bc:=u(0,t)=0,u(1,t)=0:
ic:=u(x,0)=x*(1-x):
sol:=pdsolve({pde,ic,bc},u(x,t)):

solA:=subs(sum=Sum,value(eval(eval(sol,t=2),Int=int))) assuming x>0, x<1;

plot(rhs(solA),x=0..1) assuming x>0, x<1;