# Best method to solve this discretized PDE

• I
• Daniel Sellers
In summary, the conversation discusses an attempt to numerically solve a PDE involving the function Ψ, representing a stream function on a 2D annulus grid. The conversation also mentions the use of a vertex centered discretization and an initial guess method to solve the problem. Additionally, the conversation explores the idea of treating the known terms in the equation as a group and using the Guass-Seidel method to find the solution. Finally, the conversation mentions the potential use of the analytical solution for a constant density to test the numerical solution.

#### Daniel Sellers

TL;DR Summary
I have a discretized second order PDE which I am attempting to solve numerically and having a hard time. Looking for suggestions for the best way to solve it.
I am attempting to solve the following PDE for Ψ representing a stream function on a 2D annulus grid:

(1/s)⋅(∂/∂s)[(s/ρ)(∂ψ/∂s)] + (1/s2)⋅(∂/∂Φ)[(1/ρ)(∂ψ/∂Φ)] - 2Ω + ρ(c0 + c1ψ) = 0

I have made a vertex centered discretization:

(1/sj)⋅(1/Δs2)⋅[(sj+1/2j+1/2,l){ψj+1,l - ψj,l} - (sj-1/2j-1/2,l){ψj,l - ψj-1,l}] + (1/sj2)⋅(1/ΔΦ2)⋅[(1/ρj,l+1)⋅{ψj,l+1 - ψj,l} - (1/ρj,l-1)⋅{ψj,l- ψj,l-1}] - 2Ω + ρj,l(c0 + c1Ψj,l)

This can be made to form a tridiagonal matrix by iterating through j with constant l, if we can treat ψj,l+1 and ψj,l-1 terms as known values. It's been suggested to me that we can group these two terms together (R = k1ψj,l+1 + k2ψj,l-1) and solve for them based on an 'initial guess' of all Ψj,l=c.

The resulting tridiagonal matrix (in terms of ψj,l, ψj-1,l, ψj+1,l) can then be solved as long as the boundary conditions are specified. The problem I'm running into is how to actually treat R.

Is it valid to assume R = k1ψj,l+1 + k2ψj,l-1 + 2Ω - ρj,lc0 (group all of the known terms in the equation together)? If the last two terms are included in R then the method just results in solving for the initial guess given.

Alternatively, we could base the R1 on an initial guess, solve for Ψ1,l in terms of Ψ2,l and then find R2 in terms of Ψ1,l, Ψ2,l, Ψ3,l. This let's us find Ψ2,l in terms of only Ψ3,l and so on. I believe this is similar methodology to the Guass-Seidel method of updating in place.

Neither of these methods seems to converge over many iterations. The first method, where R is based completely on a 'guess' vector with 2Ω - ρj,lc0 implicitly included, the values for psi simply decrease slowly over each iteration. The second method I've described appears to be even more unstable and quickly blows up to huge values no matter the initial guess.

Thanks for reading this far, I realize my approach may be somewhat naive as I have not developed intuition for this kind of problem solving yet. How would you solve this PDE numerically?

#### Attachments

• Vertex_Center.png
4.9 KB · Views: 403
Delta2
I don't understand your problem. You apparently have a function ψ which is a function of s and φ. But then what is ρ? Is ρ a given function? Is ρ the radial coordinate? Is there some relation between s and ρ? Please explain further.

Sorry I did neglect to identify my variables in that long post.
ρ is mass density, assumed to be a known function of s and Φ.
The work I've done so far trying to solve this numerically has been in python.

I have some observations about your problem. First, this problem has an analytic solution. Consider the mass flux vector ##\rho \vec V##; for low speed flows ##\rho## can be considered constant and one makes the substitution in the PDE$$\psi (r,\theta) = \frac \psi \rho$$ Write the PDE $$\frac 1 r \frac {\partial} {\partial r} (r\frac {\partial \psi} {\partial r} ) +\frac 1 {r^2} \frac {\partial^2 \psi} {\partial \theta^2}= a_0 + b\psi$$ where $$a_0=2\Omega-\rho c_0\\ b= \rho c_1$$Make the substitutions(not differentiation)$$\psi'=\psi - a_0\\ \psi'=R(r) \Theta(\theta)$$now use "separation of variables" and using primes to indicate differentiation, we have$$r^2 \frac {R' '} {R} + r\frac {R'} {R}-r^2b =-\frac {\Theta' '} {\Theta}$$The dependent variable of the l.h.s. is independent of the r.h.s. so both l.h.s. and r.h.s. must be equal to a constant, k^2. We get two second order differential equations$$r^2 R' ' + rR' +(r^2b-k^2)R=0\\\Theta' ' +k^2\Theta=0$$The first equation is known as Bessel's parametric equation with the solution, ##R(r)=J_k(\sqrt b r) + Y_k(\sqrt b r)## . The second equation has the solution ##\Theta (\theta)=A\cos (k\theta) + B\sin (k\theta)## where A and B are constants to be determined by boundary values. We thus have for the stream function,$$\psi (r,\theta)=\rho(J_k(\sqrt b r) + Y_k(\sqrt b r))(A\cos (k\theta) + B\sin (k\theta)) + \rho a_0$$If, however, this problem is an exercise in numerical analysis, I suggest that you first separate variables as I have done and when laying out your grid points, using the numerical recipes for first and second differentiation, you observe that ##\Delta \theta## is an arc length. The problem requires you to invert two tri-diagonal matrices.

Daniel Sellers
This is a numerical problem solving exercise. You're suggesting solve the two second order ODEs for R(r) and Θ(θ) by discretizing them and generating a tridiagonal matrix?

I'm not sure I can justify defining ρ as constant in the first step though. We are assuming equilibrium (dψ/dt) = 0 but the professor who put me onto this problem explicitly said that we have to treat ρ as dependent on the spatial coordinates.

As an update if anyone is still interested, I constructed a guass-seidel algorithm which converges for most of the examples I've tried.
One of my next steps will be to test the solution for a constant density using the analytical solution, so thanks for helping with that!

Fred Wright said:
I have some observations about your problem. First, this problem has an analytic solution. Consider the mass flux vector ##\rho \vec V##; for low speed flows ##\rho## can be considered constant and one makes the substitution in the PDE$$\psi (r,\theta) = \frac \psi \rho$$ Write the PDE $$\frac 1 r \frac {\partial} {\partial r} (r\frac {\partial \psi} {\partial r} ) +\frac 1 {r^2} \frac {\partial^2 \psi} {\partial \theta^2}= a_0 + b\psi$$ where $$a_0=2\Omega-\rho c_0\\ b= \rho c_1$$Make the substitutions(not differentiation)$$\psi'=\psi - a_0\\ \psi'=R(r) \Theta(\theta)$$now use "separation of variables" and using primes to indicate differentiation, we have$$r^2 \frac {R' '} {R} + r\frac {R'} {R}-r^2b =-\frac {\Theta' '} {\Theta}$$The dependent variable of the l.h.s. is independent of the r.h.s. so both l.h.s. and r.h.s. must be equal to a constant, k^2. We get two second order differential equations$$r^2 R' ' + rR' +(r^2b-k^2)R=0\\\Theta' ' +k^2\Theta=0$$The first equation is known as Bessel's parametric equation with the solution, ##R(r)=J_k(\sqrt b r) + Y_k(\sqrt b r)## . The second equation has the solution ##\Theta (\theta)=A\cos (k\theta) + B\sin (k\theta)## where A and B are constants to be determined by boundary values. We thus have for the stream function,$$\psi (r,\theta)=\rho(J_k(\sqrt b r) + Y_k(\sqrt b r))(A\cos (k\theta) + B\sin (k\theta)) + \rho a_0$$If, however, this problem is an exercise in numerical analysis, I suggest that you first separate variables as I have done and when laying out your grid points, using the numerical recipes for first and second differentiation, you observe that ##\Delta \theta## is an arc length. The problem requires you to invert two tri-diagonal matrices.
A couple of things I'd like to clarify as I'm having trouble with this. Shouldn't b = ρ2c1 per the first substitution?

For the second set of substitutions, if we replace Ψ in the eq by Ψ' I can't get a0 to vanish algebraically as you have.

Sorry if this is obvious, I've only worked with separation of variables in a couple of limited contexts.