6
$\begingroup$

I wanna solve the following heat transfer PDE using Mathematica.

$\qquad u_{xx}=u_{t}$

with following conditions:

$\qquad \begin{cases}u(x,0)=sin(x) &0<x<\pi &,t>0\\u_{x}(0,t)=1\\u_{x}(\pi,t)=-1\end{cases}$

I used the following codes:

heqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == Sin[x];
bc = {Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1};
sol = First @ DSolve[{heqn, ic, bc}, u[x, t], {x, t}]
sol = u[x, t] /. sol /. {K[1] -> n, Infinity -> 20}

But Mathematica says:

ReplaceAll::rmix: Elements of {(u^(0,1))[x,t]==(u^(2,0))[x,t],u[x,0]==Sin[x],{(u^(1,0))[0,t]==1,(u^(1,0))[[Pi],t]==-1}} are a mixture of lists and nonlists.

Where is the problem?

$\endgroup$
  • $\begingroup$ use bc = Sequence[Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1] to fix the error indicated by the message. $\endgroup$ – kglr Jan 17 at 19:33
  • $\begingroup$ @Nasser sorry i didn't understand what do you mean. Dsolve just gives the PDE in your image. Where is the solution? $\endgroup$ – Moreza7 Jan 17 at 19:36
  • $\begingroup$ @kglr No difference. $\endgroup$ – Moreza7 Jan 17 at 19:40
  • $\begingroup$ Kevin, do you still get the same error message? $\endgroup$ – kglr Jan 17 at 19:45
  • $\begingroup$ @Nasser How about using numerical method?If it does not have a solution, so why when I use NDSolve and plot the results, it draws a plot? Sorry I'm not pro at this topic. $\endgroup$ – Moreza7 Jan 17 at 19:59
9
$\begingroup$

I solved this by hand to confirm Maple solution.

Since the boundary conditions are not homogeneous, we can't use separation of variables. Let the solution be

$$ u=v\left( x,t\right) +r\left( x\right) $$ Where $v\left( x,t\right) $ is the solution to $v_{t}=v_{xx}$ and homogenous B.C. $v_{x}\left( 0,t\right) =0,v_{x}\left( \pi,t\right) =0$ and $r\left( x\right) $ is any reference solution which only needs to satisfy the nonhomogeneous boundary conditions: $r^{\prime}\left( 0\right) =1,r^{\prime }\left( \pi\right) =-1$. By guessing, let $r\left( x\right) =Ax+Bx^{2}$. Let see if this satisfies the boundary conditions. $r^{\prime}=A+2Bx$. At $x=0$ this implies $1=A$. Hence $r=x+Bx^{2}$. Now $r^{\prime}=1+2Bx$. At $x=\pi$ this gives $-1=1+2B\pi$ or $B=-\frac{1}{\pi}$. Therefore $$ r\left( x\right) =x-\frac{1}{\pi}x^{2} $$ Substituting $u=v\left( x,t\right) +r\left( x\right) $ into the PDE\ $u_{t}=u_{xx}$ and noting that $r^{\prime\prime}\left( x\right) =-\frac{2}{\pi}$ gives \begin{equation} v_{t}=v_{xx}-\frac{2}{\pi}\tag{1} \end{equation} PDE (1) is now solved using eigenfunction expansion. We need to find eigenfunctions and eigenvalues of $v_{t}=v_{xx}$ with $v_{x}\left( 0,t\right) =0,v_{x}\left( \pi,t\right) =0$. This is known PDE and have eigenfunctions and eigenvalues as follows. For zero eigenvalue, the eigenfunction is an arbitrary constant. Say $\beta$. let $\beta=1$ since scale is not important. $$ \Phi_{0}\left( x\right) =1 $$ And for $n=1,2,3,\cdots$ \begin{align*} \Phi_{n}\left( x\right) & =\cos\left( \sqrt{\lambda_{n}}x\right) \\ & =\cos\left( nx\right) \end{align*} with eigenvalues $\lambda_{n}=n^{2}$ for $n=1,2,3,\cdots$. Now we can eigenfunction expansion and assume the solution to (1)\ is \begin{equation} v\left( x,t\right) =\sum_{n=0}^{\infty}A_{n}\left( t\right) \Phi _{n}\left( x\right) \tag{2} \end{equation} Plugging this into the PDE (1) gives $$ \sum_{n=0}^{\infty}A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) =\sum_{n=0}^{\infty}A_{n}\left( t\right) \Phi_{n}^{\prime\prime}\left( x\right) -\frac{2}{\pi} $$ But $\Phi_{n}^{\prime\prime}\left( x\right) =-\lambda_{n}\Phi_{n}\left( x\right) $ and the above simplifies to $$ \sum_{n=0}^{\infty}A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) =-\sum_{n=0}^{\infty}A_{n}\left( t\right) \lambda_{n}\Phi_{n}\left( x\right) -\frac{2}{\pi} $$ Since eigenfunctions are complete, we can expand $\frac{2}{\pi}$ using them and the above becomes \begin{align} \sum_{n=0}^{\infty}A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) & =-\sum_{n=0}^{\infty}A_{n}\left( t\right) \lambda_{n}\Phi_{n}\left( x\right) -\sum_{n=0}^{\infty}C_{n}\Phi_{n}\left( x\right) \nonumber\\ A_{n}^{\prime}\left( t\right) \Phi_{n}\left( x\right) +A_{n}\left( t\right) \lambda_{n}\Phi_{n}\left( x\right) & =-C_{n}\Phi_{n}\left( x\right) \nonumber\\ A_{n}^{\prime}\left( t\right) +A_{n}\left( t\right) \lambda_{n} & =-C_{n}\tag{3} \end{align} To find $C_{n}$ $$ \sum_{n=0}^{\infty}C_{n}\Phi_{n}\left( x\right) =\frac{2}{\pi} $$ For $n=0$ $$ C_{0}\Phi_{0}\left( x\right) =\frac{2}{\pi} $$ But $\Phi_{0}\left( x\right) =1$, hence $$ C_{0}=\frac{2}{\pi} $$ All other $C_{m}\,\ $\ for $m>0$ are zero. Hence (3) becomes, for $n=0$ (since $\lambda_{0}=0$) \begin{align*} A_{0}^{\prime}\left( t\right) & =-\frac{2}{\pi}\\ A_{0}\left( t\right) & =-\frac{2}{\pi}t+B_{0} \end{align*} Where $B_{0}$ is integration constant. For $n>0$ (3) becomes $$ A_{n}^{\prime}\left( t\right) +A_{n}\left( t\right) n^{2}=0 $$ This has the solution $$ A_{n}\left( t\right) =B_{n}e^{-n^{2}t} $$ Where $B_{n}$ is constant of integration. Hence from (2) \begin{align*} v\left( x,t\right) & =\sum_{n=0}^{\infty}A_{n}\left( t\right) \Phi _{n}\left( x\right) \\ & =A_{0}\left( t\right) +\sum_{n=1}^{\infty}A_{n}\left( t\right) \Phi _{n}\left( x\right) \\ & =-\frac{2}{\pi}t+B_{0}+\sum_{n=1}^{\infty}B_{n}e^{-n^{2}t}\cos\left( nx\right) \end{align*} Since $u=v\left( x,t\right) +r\left( x\right) $ then the solution becomes \begin{equation} u\left( x,t\right) =\left( x-\frac{1}{\pi}x^{2}\right) -\frac{2}{\pi }t+B_{0}+\sum_{n=1}^{\infty}B_{n}e^{-n^{2}t}\cos\left( nx\right) \tag{4} \end{equation} At $t=0$ \begin{equation} \sin\left( x\right) =\left( x-\frac{1}{\pi}x^{2}\right) +B_{0}+\sum _{n=1}^{\infty}B_{n}\cos\left( nx\right) \tag{5} \end{equation} case $n=0$ $$ \int_{0}^{\pi}\sin\left( x\right) \cos\left( \sqrt{\lambda_{0}}x\right) dx=\int_{0}^{\pi}\left( x-\frac{1}{\pi}x^{2}\right) \cos\left( \sqrt{\lambda_{0}}x\right) dx+\int_{0}^{\pi}B_{0}\cos\left( \sqrt {\lambda_{0}}x\right) dx $$ But $\lambda_{0}=0$ hence \begin{align*} \int_{0}^{\pi}\sin\left( x\right) dx & =\int_{0}^{\pi}\left( x-\frac {1}{\pi}x^{2}\right) dx+\int_{0}^{\pi}B_{0}dx\\ 2 & =\frac{\pi^{2}}{6}+B_{0}\pi\\ B_{0} & =\frac{2}{\pi}-\frac{\pi}{6} \end{align*} For $n>0$, Multiplying both sides of (5) by $\cos\left( mx\right) $ and integrating $$ \int_{0}^{\pi}\sin\left( x\right) \cos\left( mx\right) dx=\int_{0}^{\pi }\left( x-\frac{1}{\pi}x^{2}\right) \cos\left( mx\right) dx+\sum _{n=1}^{\infty}B_{n}\int_{0}^{\pi}\cos\left( nx\right) \cos\left( mx\right) dx $$ For $m=1$ \begin{align*} 0 & =0+B_{1}\frac{\pi}{2}\\ B_{1} & =0 \end{align*} For $m>1$ \begin{align*} -\frac{1+\left( -1\right) ^{m}}{m^{2}\left( -1+m^{2}\right) } & =\frac{\pi}{2}B_{m}\\ B_{m} & =\frac{-2}{\pi}\left( \frac{1}{m^{2}}\frac{\left( -1\right) ^{m}+1}{m^{2}-1}\right) \end{align*} Hence solution (4) becomes \begin{align*} u\left( x,t\right) & =\left( x-\frac{1}{\pi}x^{2}\right) -\frac{2}{\pi }t-\frac{\pi}{6}+\frac{2}{\pi}+\sum_{n=1}^{\infty}B_{n}e^{-n^{2}t}\cos\left( nx\right) \\ u\left( x,t\right) & =\left( x-\frac{1}{\pi}x^{2}\right) -\frac{2}{\pi }t-\frac{\pi}{6}+\frac{2}{\pi}+\sum_{n=2}^{\infty}\frac{-2}{\pi}\left( \frac{1}{n^{2}}\frac{\left( -1\right) ^{n}+1}{n^{2}-1}\right) e^{-n^{2} t}\cos\left( nx\right) \end{align*}

Maple solution verified:

heqn := diff(u(x, t), t) = diff(u(x, t), x$2):
ic := u(x, 0) = sin(x):
bc := eval(diff(u(x,t),x),x=0)=1,  eval( diff(u(x,t),x),x=Pi)=-1:
sol := pdsolve({heqn, ic, bc}, u(x, t))

enter image description here

Verification of hand solution (used 10 terms in sum, good enough):

mysol=(x-1/Pi x^2)-2/Pi t-Pi/6+2/Pi-2/Pi 
           Sum[((-1)^n+1)/(n^2(n^2-1)) Cos[n x]Exp[-n^2 t],{n,2,10}];
D[mysol,x]/.t->0/.x->0
(*1*)
D[mysol,x]/.t->0/.x->Pi
(*-1*)
Plot[{mysol/.t->0},{x,0,Pi},AxesOrigin->{0,0}]

Mathematica graphics

Animation of hand solution

Manipulate[
 Plot[Evaluate[mysol /. t -> time], {x, 0, Pi}, 
  PlotRange -> {{0, Pi}, {-0.5, 1}}],
 {{time, 0, "time"}, 0, 1, .1},
 TrackedSymbols :> {time},
 Initialization :> {mysol = (x - 1/Pi x^2) - 2/Pi t - Pi/6 + 2/Pi - 
     2/Pi Sum[((-1)^n + 1)/(n^2 (n^2 - 1)) Cos[n x] Exp[-n^2 t], {n, 2, 
         10}]}
 ]

enter image description here

ref (1): boundary_values

$\endgroup$
  • 1
    $\begingroup$ 1++,Great answer . $\endgroup$ – Mariusz Iwaniuk Jan 18 at 5:48
  • 1
    $\begingroup$ Changing the coefficient of t is not the same as adding any number to the solution. And indeed, changing that coefficient does not effect the bc's or ic's, but it breaks the pde. $\endgroup$ – Bill Watts Jan 18 at 7:36
  • 1
    $\begingroup$ The premise that only having neumann boundary conditions makes the solution non-unique does not apply to the heat equation. There are numerous examples all over the web. One is found at ramanujan.math.trinity.edu/rdaileda/teach/s12/m3357/lectures/…. For some pde's the premise is correct, but different pde's have different bc requirements. By the way, separation of variables works very well in this case, but you do have to use plus separation as well as times separation to get all the pieces. $\endgroup$ – Bill Watts Jan 18 at 19:42
  • $\begingroup$ The B1 must be zero. How it is 2/Pi?! $\endgroup$ – Moreza7 Jan 21 at 20:30
  • $\begingroup$ @Nasser Thank You. But using Mathematica for Fourier series coefficient: FourierCosCoefficient[Sin[x] - (x - ((x^2)/Pi)), x, 0] gets: 4/Pi - Pi/3 which is twice as yours: 2/Pi-P/6 why? $\endgroup$ – Moreza7 Jan 21 at 20:50
6
$\begingroup$
 MAPLESOL = (x - 1/Pi x^2) - 2/Pi t - Pi/6 + 2/Pi - 2/Pi Inactivate[
 Sum[((-1)^n + 1)/(n^2 (n^2 - 1)) Cos[n x] Exp[-n^2 t], {n, 2, Infinity}]];

 Plot3D[MAPLESOL /. Infinity -> 20 // Activate, {x, 0, Pi}, {t, 0, 1}]

Check equation,initial and boundary conditions:

 (D[MAPLESOL /. Infinity -> 20, t] == D[MAPLESOL /. Infinity -> 20, {x, 2}]) // Activate
 (*True*)

 D[MAPLESOL /. Infinity -> 20 // Activate, x] /. t -> 0 /. x -> 0
 (* 1 *)
 D[MAPLESOL /. Infinity -> 20 // Activate, x] /. t -> 0 /. x -> Pi
 (* -1 *)

 Plot[{Sin[x], Evaluate[(MAPLESOL /. Infinity -> 20 // Activate) /. t -> 0]}, {x, 0, Pi}, PlotStyle -> {Red, {Dashed, Black}}]

 Plot3D[D[MAPLESOL /. Infinity -> 20 // Activate, x] // Evaluate, {x, 0, Pi}, {t, 0, 1}]

It's seems symbolic solution by Maple is OK.

enter image description here

$\endgroup$
  • $\begingroup$ good answer. Yes, this seems like a valid solution. But it is not unique solution as adding any number to the solution is also a solution. D[(10 + MAPLESOL) /. Infinity -> 20 // Activate, x] /. t -> 0 /. x -> 0 gives 1 as well. May be Mathematica did not solve it because of this? I do not know. Maple 2018 generally can solve more PDE's than Mathematica 11.3 $\endgroup$ – Nasser Jan 18 at 0:57
  • $\begingroup$ Adding a constant to the stated solution will break the ic u[x,0] = Sin[x]. I think that the given ic is enough to make the solution unique.When you work the solution with separation of variables there are no extra constants lying around. $\endgroup$ – Bill Watts Jan 18 at 1:10
  • $\begingroup$ @Nasser and Mariusz …This solution is incorrect, just compare it to the numeric solution. $\endgroup$ – xzczd Jan 18 at 10:30
5
$\begingroup$

This question is strongly related to, if not a duplicate of this question. It can be solved with the help of finite Fourier cosine transform and its inversion:

heqn = D[u[x, t], t] == D[u[x, t], {x, 2}];
ic = u[x, 0] == Sin[x];
bc = {Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1};

help[index_] := 
 Module[{tset = 
    finiteFourierCosTransform[{heqn, ic}, {x, 0, Pi}, index] /. Rule @@@ bc /. 
     HoldPattern@finiteFourierCosTransform[f_, __] :> f}, 
  tsol = DSolve[tset, u[x, t], t][[1, 1, -1]]]

tsolgeneral = help[n] // Simplify

tsolzero = help[0]

tsolfunc[n_] = Piecewise[{{tsolgeneral, n != 0}}, tsolzero]

sol = inverseFiniteFourierCosTransform[tsolfunc[n], n, {x, 0, Pi}]

sol = sol /. HoldForm@Sum[expr_, {n, C}] :> sum[Simplify[expr, n > 1], {n, 2, C}] /. 
 sum -> (HoldForm@Sum@## &)

Mathematica graphics

Notice finite Fourier cosine transform for $n=0$ has been calculated separately because currently finiteFourierCosTransform, which is built on Integrate, isn't strong enough to obtain the complete solution for general n.

Finally let's compare the solution to the numeric one:

solapprox = Compile[{t, x}, #] &[sol /. C -> 10 // ReleaseHold];

solnum = NDSolveValue[{heqn, ic, bc}, u, {t, 0, 1}, {x, 0, Pi}]

Manipulate[Plot[{solnum[x, t], solapprox[t, x]}, {x, 0, Pi}], {t, 0, 1}]

enter image description here


Addendum

The solution found above looks different from the one found by method of separation of variables shown in other answers, but it's possible to prove they're identical i.e.

sol= -((2 (-1 + t))/Pi) + (2 
    HoldForm[Sum[-(((1 + (-1)^n) (1 + E^(n^2 t) (-1 + n^2)) 
           Cos[n x])/(E^(n^2 t) (n^2 (-1 + n^2)))), {n, 2, C}]])/Pi

is equivalent to

mysol=(x-1/Pi x^2)-2/Pi t-Pi/6+2/Pi-
   2/Pi Inactivate[Sum[((-1)^n+1)/(n^2(n^2-1)) Cos[n x]Exp[-n^2 t],{n,2,Infinity}]];

We first find the difference between the summand and calculate the sum:

summandseparate = (((-1)^n + 1) Cos[n x] Exp[-n^2 t])/(n^2 (n^2 - 1));

summandtransform = ((1 + (-1)^n) E^(-n^2 t) (1 + E^(n^2 t) (-1 + n^2)) Cos[n x])/(
  n^2 (-1 + n^2));

diff = Sum[summandseparate - summandtransform, {n, 2, Infinity}] // FullSimplify
(* 1/4 (-PolyLog[2, E^(-2 I x)] - PolyLog[2, E^(2 I x)]) *)

Remark

At least in v11.2, you need to modify definition of diff to

diff = Sum[summandseparate - summandtransform // Simplify // Evaluate, {n, 2, Infinity}]

or the calculation will be extremely slow.

Then the problem boils down to proving the following identity:

eq = (x - x^2/π) - (2 t)/π - π/6 + 2/π - 2/Pi diff2 == -((
    2 (-1 + t))/π);

-4 diff == -4 diff2 /. First@Solve[eq, diff2]
(* PolyLog[2, E^(-2 I x)] + PolyLog[2, E^(2 I x)] == 1/3 (π^2 - 6 π x + 6 x^2) *)

which seems to be beyond the reach of Mathematica, so I asked it in math.SE and get the solution by hand in 6 minutes. Please check this post for more details.

$\endgroup$
  • $\begingroup$ -4 Series[diff, {x, 0, 10}] // FullSimplify[#, 0 < x < \[Pi]] & gives you what you want. I guess your tails are harder to get rid of because your series converges much more slowly than the others because your Exp[-n^2 t] only damps out part of the sum while the sum damps out very rapidly in the separate variable case. I still don't understand why I can't prove your solution satisfies the pde. Probably need to know a closed form for the whole sum to do that. $\endgroup$ – Bill Watts Jan 20 at 0:53
  • 1
    $\begingroup$ @Bill Yeah Mathematica can verify the identity in experimental mathematics style. We can prove the solution satisfy the PDE in the same manner as shown in Addendum, the key point is to sum part of the summand, the following code takes about 100 seconds: terms = List @@ sol /. HoldForm@Sum[a_, _] :> Sequence@a; sum = Sum[#, {n, 2, Infinity}] &; {Simplify@D[terms, t], {#, sum /@ Expand@#2} & @@ D[terms, x, x] // Simplify} // AbsoluteTiming $\endgroup$ – xzczd Jan 20 at 6:32
  • 1
    $\begingroup$ See my update. We don't need the math guys for this. $\endgroup$ – Bill Watts Jan 20 at 22:23
2
$\begingroup$

And yet another way.

pde = D[u[x, t], t] == D[u[x, t], {x, 2}]

ic = u[x, 0] == Sin[x]

bc = {Derivative[1, 0][u][0, t] == 1, Derivative[1, 0][u][Pi, t] == -1}

Separate variables, first with times separation.

u[x_, t_] = X[x] T[t]

pde/u[x, t] // Expand
T'[t]/T[t] == X''[x]/X[x]

Each side must be equal to a constant. First try 0.

t0eq = T'[t]/T[t] == 0;

DSolve[t0eq, T[t], t] // Flatten
(*{T[t] -> C[1]}*)

t0 = 1

since we will combine with constants from the x equation.

x0eq = X''[x]/X[x] == 0

DSolve[x0eq, X[x], x] // Flatten
(*{X[x] -> C[2] x + C[1]}*)

x0 = X[x] /. % /. {C[1] -> c1, C[2] -> c2}

Now use a negative constant.

t1eq = T'[t]/T[t] == -α^2

DSolve[t1eq, T[t], t] // Flatten
(*{T[t] -> C[1] E^(α^2 (-t))}*)

t1 = T[t] /. % /. C[1] -> 1

x1eq = X''[x]/X[x] == -α^2

DSolve[x1eq, X[x], x] // Flatten
(*{X[x] -> C[2] Sin[α x] + C[1] Cos[α x]}*)

x1 = X[x] /. % /. {C[1] -> c3, C[2] -> c4}

We don't have all the necessary pieces yet, so do plus separation.

u[x_, t_] = X[x] + T[t]

pde
(*T'[t] == X''[x]*)

Again, each side must be equal to a constant. Call it δ.

xpeq = X''[x] == δ

DSolve[xpeq, X[x], x] // Flatten
(*{X[x] -> C[2] x + C[1] + (δ x^2)/2}*)

xp = X[x] /. % /. {C[1] -> c5, C[2] -> c6}

tpeq = T'[t] == δ

DSolve[tpeq, T[t], t] // Flatten
(*{T[t] -> C[1] + δ t}*)

tp = T[t] /. % /. C[1] -> 0

Put all the pieces together.

Clear[c1, c2, c3, c4, c5, c6, α, δ]

u[x_, t_] = x0 t0 + x1 t1 + xp + tp
(*c1 + c2 x + E^(α^2 (-t)) (c3 Cos[α x] + c4 Sin[α x]) + c5 + 
 c6 x + δ t + (δ x^2)/2*)

c1 and c5 can be combined, c2 and c6 can be combined.

c5 = 0;
c6 = 0;

Look at the first bc

bc[[1]]
(*c2 + α c4 E^(α^2 (-t)) == 1*)

from which we can say

c2 = 1;
c4 = 0;

and the next bc.

bc[[2]]
-(α c3 Sin[π α] E^(α^2 (-t))) + π δ + 1 == -1

Make the Sin 0

Solve[α π == n π, α] // Flatten

from which

α = n;
$Assumptions = n ∈ Integers && n >= 0

Use the rest of the equation to solve for δ

c26eq2 = π δ + 1 == -1

δ = δ /. Solve[π δ + 1 == -1, δ][[1]]

So far we have

u[x, t]
c1 + c3 E^(-n^2 t) Cos[n x] - (2 t)/π - x^2/π + x

Now for the ic.

ic
(*c1 + c3 Cos[n x] - x^2/π + x == Sin[x]*)

We are going to have an infinite series for n and the n=0 will be a constant term. Since we don't need two of them we can throw out c1.

c1 = 0

ic - (-(x^2/π) + x)
(*c3 Cos[n x] == x^2/π - x + Sin[x]*)

To solve for c3, multiply by Cos[n x] and integrate.

c3*Integrate[Cos[n*x]^2, {x, 0, Pi}] == (1/Pi)*Integrate[x^2*Cos[n*x], {x, 0, Pi}] - 
   Integrate[x*Cos[n*x], {x, 0, Pi}] + Integrate[Sin[x]*Cos[n*x], {x, 0, Pi}]

c3 = c3 /. Solve[%, c3][[1]] // Simplify
(*-((2 ((-1)^n + 1))/(\[Pi] n^2 (n^2 - 1)))*)

The constant term for n = 0 needs to be done separately.

c30*Integrate[1, {x, 0, Pi}] == (1/Pi)*Integrate[x^2*1, {x, 0, Pi}] - 
   Integrate[x*1, {x, 0, Pi}] + Integrate[Sin[x]*1, {x, 0, Pi}]

c30 = c30 /. Solve[%, c30][[1]] // Expand
2/π-π/6

Without the sum on n we have

u[x_, t_] = u[x, t] + c30
-((2 ((-1)^n + 1) E^(-n^2 t) Cos[n x])/(π n^2 (n^2 - 1))) - (
 2 t)/π - x^2/π + x - π/6 + 2/π

Check against the pde.

pde
(*True*)

We can see that odd n the terms are zero, so change to.

-((2 ((-1)^n + 1) E^(-n^2 t) Cos[n x])/(π n^2 (n^2 - 1))) /. 
  n -> 2 n // Simplify
(*-((E^(-4 n^2 t) Cos[2 n x])/(π n^2 (4 n^2 - 1)))*)

We have an infinite series of n starting from 1, but lets make a finite series.

u[x_, t_, mm_] := 2/Pi - Pi/6 - (2*t)/Pi - x^2/Pi + x - 
   (1/Pi)*Sum[Cos[2*n*x]/(E^(4*n^2*t)*(n^2*(4*n^2 - 1))), {n, 1, mm}]

Check against xzczd's numerical solution.

gifs = Table[
   Plot[{u[x, t, 20], solnum[x, t]}, {x, 0, π}, 
    PlotRange -> {-1, 1}], {t, 0, 1, .01}];
ListAnimate[%]

This solution matches Maple's and Nasser's solution and the numerical solution. I could not get xzczd's analytic solution to satisfy the pde, which is I think is the cause for the little tails at the end points.

enter image description here

Update on showing the equivalence in the separate variable solution and xzczd's fourier transform solution. The separate variable solution simplified by letting the original n go to 2n to eliminate the odd n zero's:

usep[x_, t_] := -((2*(t - 1))/Pi) - Pi/6 - x^2/Pi + x - 
   (1/Pi)*Sum[Cos[2*n*x]/(E^(4*n^2*t)*(n^2*(4*n^2 - 1))), {n, 1, Infinity}]

Rewrite as

usep[x, t] := 2/Pi - (2*t)/Pi - (1/Pi)*Sum[Cos[2*n*x]/(E^(4*n^2*t)*(n^2*(4*n^2 - 1))), 
     {n, 1, Infinity}] + Sum[An*Cos[2*n*x], {n, 1, Infinity}]

where we find a series such that for 0 <= x <= π

x - x^2/Pi - Pi/6 == Sum[An*Cos[2*n*x], {n, 1, Infinity}];

Find An by multiplying by Cos[2 n x] and integrating.

Integrate[(x - x^2/Pi - Pi/6)*Cos[2*n*x], {x, 0, Pi}] == 
  An*Integrate[Cos[2*n*x]^2, {x, 0, Pi}]

An = An /. Solve[%, An][[1]] 
(*-(1/(π n^2))*)

Add to the existing sum

-(1/π) (E^(-4 n^2 t) Cos[2 n x])/(n^2 (4 n^2 - 1)) + 
  An Cos[2 n x] // Simplify
(*((E^(-4 n^2 t)/(1 - 4 n^2) - 1) Cos[2 n x])/(π n^2)*)

And we end up with the same solution as the fourier series solution.

$\endgroup$
  • $\begingroup$ There're little tails at the end points in the graphic of solapprox because I've truncated the summation at $n=10$ (C -> 10). If e.g. C -> 60 is used, the tails will be almost invisible. And, it's possible to prove my solution is identical to the one found by method of separation of variables, see my update. $\endgroup$ – xzczd Jan 19 at 12:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.