Modeling a Partial Differential Equation

I'm in the midst of solving the given problem

Given the setup I tried writing up the following:

ClearAll["Global*"]
β = 0.05;
sol = NDSolveValue[{
D[Subscript[U, z][η, τ], τ] == -Sin[τ] + β (1/η D[η*Subscript[U, z][η, τ], η, η]),
Subscript[U, z][η, 0] == 0, Subscript[U, z][1, τ] == 0,},
Subscript[U, z], {τ, 0, 10}]


How can I model the second boundary condition accordingly, and what should I do next to model the solution on a 3D plot?

• Have a look at Menu/Help/WolframDocumentation/NeumannValue and there go to Applications/Time-Dependent Problems. There you will find examples which are no exactly equal to yours (that is, not radial), but very close and with comparable type of Neumann and initial conditions. Nov 19 '19 at 8:48
• The book does not reference that function so I will not use that Nov 19 '19 at 21:22

A couple things. Since $$\eta$$ is in the denominator of your pde, your solution will have a problem when $$\eta$$ is 0. Luckily since zero is one of the end points we can pick a number close to zero without changing the answer much. A numerical solution is an approximation after all. Also your NDSolveValue needs the limits of $$\eta$$ as well as $$\tau$$. Fixing that and choosing a very small number for zero, we have:

ClearAll["Global*"]
β = 0.05;
zero = 10^-12;

sol = NDSolveValue[{D[U[η, τ], τ] == -Sin[τ] + β*(1/η)*D[η*D[U[η, τ], η], η],
U[η, 0] == 0, U[1, τ] == 0, Derivative[1, 0][U][zero, τ] == 0}, U, {η, zero, 1},
{τ, 0, 10}]

tp = Table[
Plot[Evaluate[sol[η, τ]], {η, zero, 1},
PlotRange -> {-2, 1}], {τ, 0, 10, .1}];
ListAnimate[tp]


In this particular case we can include $$\eta = 0$$ in the solution because the derivative wrt to $$\eta$$ is zero there. Expanding the pde we have

Derivative[0, 1][U][η, τ] == (β*Derivative[1, 0][U][η, τ])/η +
β*Derivative[2, 0][U][η, τ] - Sin[τ]


the term with $$\eta$$ in the denominator is $$0/0$$ at $$\eta = 0$$, so we can use L'Hospital's rule. Rewriting the pde with Piecewise functions to separate the zero case from the general case and solving we get:

Clear["Global*"]

pde = Derivative[0, 1][U][η, τ] ==
Piecewise[{{β*Derivative[2, 0][U][η, τ], η == 0},
{(β*Derivative[1, 0][U][η, τ])/η, True}}] + β*Derivative[2, 0][U][η, τ] -
Sin[τ]

β = .05

sol = NDSolveValue[{pde, U[η, 0] == 0, U[1, τ] == 0,
Derivative[1, 0][U][0, τ] == 0},
U, {η, 0, 1}, {τ, 0, 10}];


The resulting plot looks the same as before, but we can now start the solution at $$\eta = 0$$.

• What is the purpose of the code: β*Derivative[2, 0][U][η, τ]  Nov 22 '19 at 0:21
• It comes from L'Hospital's rule. When you have f[x]/g[x] = 0/0 at x = 0. You can often find the limit as x->0 by finding f'[x]/g'[x]. So when you have in this case an expression like u'[x]/x, find the limit at 0 by taking the numerator and denominator derivatives and get u''[x]/1`. Nov 22 '19 at 1:24