# 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. – Alexei Boulbitch Nov 19 '19 at 8:48
• The book does not reference that function so I will not use that – TexMexDex 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][η, τ]  – TexMexDex 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`. – Bill Watts Nov 22 '19 at 1:24