Stokes' second problem [closed]

I'm trying to solve this equation

sol = u[t, x] /. NDSolve[{D[u[t, x], t] == D[u[t, x], x, x], u[0, x] == 0,
u[t, 0] == sin[t], u[t, 5] == 0}, u, {t, 0, 10}, {x, 0, 5}]


however, everytime, I use NDsolve, it gives me the below error

ReplaceAll::reps: {NDSolve[{(u^(1,0))[t,x]==(u^(0,2))[t,x],u[0,x]==0,u[t,0]==sin[t],u[t,5]==0},u,{t,0,10},{x,0,5}]} is neither a list of replacement rules nor a valid dispatch table, and so cannot be used for replacing.

• sin is spelled Sin in Mathematica. Fix that and your problem goes away
– Bill
Jun 11, 2017 at 21:55

The PDE along with the initial and boundary conditions constitute the famous Stokes second problem. Here a comparison has been made among Maple, Mathematica and COMSOL solving the same model.

The main issue you are facing is the syntax for $$sin$$ which @Bill has correctly pointed out that it should be Sin with a capital S not small s.

Sometime back, I have tried to solve the same thing, here is my try, which is not different than yours,

ClearAll["Global*"];
pde = {Derivative[1, 0][U][t, x] == Derivative[0, 2][U][t, x]}
ics = {U[0, x] == 0};
bcs = {U[t, 0] == Sin[t], U[t, 10] == 0};
ibcAll = Flatten[{ics, bcs}, 1];
sol = NDSolve[{pde, ibcAll}, {U}, {t, 0, 12}, {x, 0, 10}];

DensityPlot[(U /. First@sol)[t, x], {x, 0, 4}, {t, 0, 12},
PlotRange -> All, ColorFunctionScaling -> False,
ColorFunction -> "BrightBands", PlotPoints -> 200, MaxRecursion -> 5]


After some struggle, I found an analytical solution to this problem.

$$\left( \sum_{n=1}^{\infty}\left( \frac{2L^{2}\left( kn^{2}\pi^{2}\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi^{4}\right) }e^{-k\left( \frac{n^{2}\pi^{2}}{L^{2}}\right) t}-\frac{2L^{2}\left( kn^{2}\pi^{2}\cos t+L^{2}\sin t\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi ^{4}\right) }\right) \sin\left( \frac{n\pi}{L}x\right) \right) +\left( \frac{L-x}{L}\right) \sin\left( t\right)$$

Below is the derivation, and simulation for $k=1,L=5$ showing it matches the numerical solution from NDSolve

Solve

\begin{align} \frac{\partial u}{\partial t} & =k\frac{\partial^{2}u}{\partial x^{2}} \tag{1}\\ 0 & <x<L\nonumber\\ t & >0\nonumber \end{align}

Initial conditions

$$u\left( 0,x\right) =0$$

Boundary conditions

\begin{align*} u\left( 0,t\right) & =\sin\left( t\right) \\ u\left( t,L\right) & =0 \end{align*}

Let $$u=v+u_{E}\tag{2}$$ where $u_{E}\left( x,t\right)$ is steady state solution that only needs to satisfy boundary conditions and $v\left( x,t\right)$ satisfies the PDE itself but with homogenous B.C. At steady state, the PDE becomes

\begin{align*} 0 & =k\frac{d^{2}u_{E}}{dx^{2}}\\ u_{E}\left( 0\right) & =\sin\left( t\right) \\ u_{E}\left( L\right) & =0 \end{align*}

The solution is $u_{E}\left( t\right) =\left( \frac{L-x}{L}\right) \sin\left( t\right)$. Hence (2) becomes

$$u\left( x,t\right) =v\left( x,t\right) +\left( \frac{L-x}{L}\right) \sin\left( t\right)$$

Substituting the above in (1) gives

\begin{align} \frac{\partial v}{\partial t}+\left( \frac{L-x}{L}\right) \cos\left( t\right) & =k\frac{\partial^{2}v}{\partial x^{2}}\nonumber\\ \frac{\partial v}{\partial t} & =k\frac{\partial^{2}v}{\partial x^{2} }+\left( \frac{x-L}{L}\right) \cos\left( t\right) \nonumber\\ \frac{\partial v}{\partial t} & =k\frac{\partial^{2}v}{\partial x^{2} }+Q\left( x,t\right) \tag{3} \end{align}

With boundary conditions $u_{0}\left( 0,t\right) =0,u\left( L,t\right) =0$. This is now in standard form and separation of variables can be used to solve it. $$Q\left( x,t\right) =\left( \frac{x-L}{L}\right) \cos\left( t\right)$$ Now acts as a source term. The eigenfunctions are known to be $\Phi_{n}\left( x\right) =\sin\left( \sqrt{\lambda_{n}}x\right)$ where $\lambda _{n}=\left( \frac{n\pi}{L}\right) ^{2}$. Hence by eigenfunction expansion, the solution to (3) is

$$v\left( x,t\right) =\sum_{n=1}^{\infty}B_{n}\left( t\right) \Phi _{n}\left( x\right) \tag{3A}$$

Substituting this into (3) gives

$$\sum_{n=1}^{\infty}\frac{dB_{n}\left( t\right) }{dt}\Phi_{n}\left( x\right) =k\sum_{n=1}^{\infty}B_{n}\left( t\right) \Phi_{n}^{\prime\prime }\left( x\right) +Q\left( x,t\right) \tag{4}$$

Expanding $Q\left( x,t\right)$ using same basis (eigenfunctions) gives

$$Q\left( x,t\right) =\sum_{n=1}^{\infty}q_{n}\left( t\right) \Phi _{n}\left( x\right)$$

Applying orthogonality

\begin{align*} \int_{0}^{L}Q\left( x,t\right) \Phi_{m}\left( x\right) dx & =\int_{0} ^{L}\sum_{n=1}^{\infty}q_{n}\left( t\right) \Phi_{n}\left( x\right) \Phi_{m}\left( x\right) dx\\ & =\sum_{n=1}^{\infty}q_{n}\left( t\right) \int_{0}^{L}\Phi_{n}\left( x\right) \Phi_{m}\left( x\right) dx \end{align*}

But $\sum_{n=1}^{\infty}\int_{0}^{L}\Phi_{n}\left( x\right) \Phi_{m}\left( x\right) dx=\int_{0}^{L}\Phi_{m}^{2}\left( x\right) dx=\frac{L}{2}$ since $\Phi_{n}\left( x\right) =\sin\left( \frac{n\pi}{L}x\right)$ and the above simplifies to

\begin{align*} \int_{0}^{L}Q\left( x,t\right) \Phi_{n}\left( x\right) dx & =\frac{L} {2}q_{n}\left( t\right) \\ q_{n}\left( t\right) & =\frac{2}{L}\int_{0}^{L}Q\left( x,t\right) \sin\left( \frac{n\pi}{L}x\right) dx \end{align*}

But $Q\left( x,t\right) =\left( \frac{x-L}{L}\right) \cos\left( t\right)$, hence

\begin{align*} q_{n}\left( t\right) & =\frac{2}{L}\int_{0}^{L}\left( \frac{x-L} {L}\right) \cos\left( t\right) \sin\left( \frac{n\pi}{L}x\right) dx\\ & =\frac{-2}{n\pi}\cos\left( t\right) \end{align*}

Therefore $Q\left( x,t\right) =\sum_{n=1}^{\infty}q_{n}\left( t\right) \Phi_{n}\left( x\right) =\sum_{n=1}^{\infty}\frac{-2}{n\pi}\cos\left( t\right) \sin\left( \frac{n\pi}{L}x\right)$ and (4) becomes

\begin{align*} \sum_{n=1}^{\infty}\frac{dB_{n}\left( t\right) }{dt}\Phi_{n}\left( x\right) & =k\sum_{n=1}^{\infty}B_{n}\left( t\right) \Phi_{n} ^{\prime\prime}\left( x\right) -\sum_{n=1}^{\infty}\frac{2}{n\pi}\cos\left( t\right) \sin\left( \frac{n\pi}{L}x\right) \\ \frac{dB_{n}\left( t\right) }{dt}\sin\left( \frac{n\pi}{L}x\right) & =kB_{n}\left( t\right) \left( -\frac{n^{2}\pi^{2}}{L^{2}}\sin\left( \frac{n\pi}{L}x\right) \right) -\frac{2}{n\pi}\cos\left( t\right) \sin\left( \frac{n\pi}{L}x\right) \\ \frac{dB_{n}\left( t\right) }{dt}+B_{n}\left( t\right) k\frac{n^{2}\pi ^{2}}{L^{2}} & =-\frac{2}{n\pi}\cos\left( t\right) \end{align*}

This is an ODE in $B_{n}\left( t\right)$ whose solution is

$$B_{n}\left( t\right) =C_{n}e^{-k\left( \frac{n^{2}\pi^{2}}{L^{2}}\right) t}-\frac{2L^{2}\left( kn^{2}\pi^{2}\cos t+L^{2}\sin t\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi^{4}\right) }$$

From (3A) $v\left( x,t\right)$ now becomes

$$v\left( x,t\right) =\sum_{n=1}^{\infty}C_{n}e^{-k\left( \frac{n^{2}\pi^{2} }{L^{2}}\right) t}\sin\left( \frac{n\pi}{L}x\right) -\frac{2L^{2}\left( kn^{2}\pi^{2}\cos t+L^{2}\sin t\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi ^{4}\right) }\sin\left( \frac{n\pi}{L}x\right) \tag{5}$$

To find $C_{n}$, from initial conditions, at $t=0$ the above becomes

$$0=\sum_{n=1}^{\infty}C_{n}\sin\left( \frac{n\pi}{L}x\right) -\frac {2L^{2}\left( kn^{2}\pi^{2}\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi ^{4}\right) }\sin\left( \frac{n\pi}{L}x\right)$$

Hence

$$C_{n}=\frac{2L^{2}\left( kn^{2}\pi^{2}\right) }{n\pi\left( L^{4}+k^{2} n^{4}\pi^{4}\right) }$$

Therefore (5) becomes

$$v\left( x,t\right) =\sum_{n=1}^{\infty}\left( \frac{2L^{2}\left( kn^{2} \pi^{2}\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi^{4}\right) }e^{-k\left( \frac{n^{2}\pi^{2}}{L^{2}}\right) t}-\frac{2L^{2}\left( kn^{2}\pi^{2}\cos t+L^{2}\sin t\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi^{4}\right) }\right) \sin\left( \frac{n\pi}{L}x\right)$$

And since $u=v+u_{E}$ then the solution is

$$u\left( x,t\right) =\left( \sum_{n=1}^{\infty}\left( \frac{2L^{2}\left( kn^{2}\pi^{2}\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi^{4}\right) }e^{-k\left( \frac{n^{2}\pi^{2}}{L^{2}}\right) t}-\frac{2L^{2}\left( kn^{2}\pi^{2}\cos t+L^{2}\sin t\right) }{n\pi\left( L^{4}+k^{2}n^{4}\pi ^{4}\right) }\right) \sin\left( \frac{n\pi}{L}x\right) \right) +\left( \frac{L-x}{L}\right) \sin\left( t\right)$$

To simulate

ClearAll[t, x, n]
k = 1; L0 = 5; max = 400;
u[x_, t_] =
Sum[(((2*L0^2*(k*n^2*Pi^2))/(n*Pi*(L0^4 + k^2*n^4*Pi^4)))*
Exp[(-k)*((n^2*Pi^2)/L0^2)*t] -
(2*L0^2*(k*n^2*Pi^2*Cos[t] + L0^2*Sin[t]))/(n*
Pi*(L0^4 + k^2*Pi^4*n^4)))*Sin[((n*Pi)/L0)*x],
{n, 1, max}] + ((L0 - x)/L0)*Sin[t];

Manipulate[Grid[{{"Analytical solution"},
{Plot[Evaluate[u[x,t]],{x,0,5},PlotRange->{{0,5},{-1.1,1.1}},
ImageSize->400]}}],
{{t,0,"t"},0,100,.01}
]


Here is the numerical solution from NDSolve using zhk code posted above to compare with

ClearAll["Global*"];
pdeset = {Derivative[1, 0][U][t, x] == Derivative[0, 2][U][t, x]}
ics = {U[0, x] == 0};
bcs = {U[t, 0] == Sin[t], U[t, 5] == 0};
ibcAll = {ics, bcs};
numericalSol = NDSolve[{pdeset, ibcAll}, U, {t, 0, 100}, {x, 0, 5}];
Manipulate[Grid[{{"Numerical solution"},
{Plot[Evaluate[U[t, x] /. numericalSol], {x, 0, 5},
PlotRange -> {{0, 5}, {-1, 1}}, ImageSize -> 400]}}],
{{t, 0, "t"}, 0, 100, .01}
]
`

• A wonderful exercise, but my choice is to first let Mathematica work. With Mathematica 11.1.1 one can choose between the numerical and analytical solution.
– user36273
Jun 12, 2017 at 9:34