Best method to solve this discretized PDE

[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/s_j)⋅(1/Δs²)⋅[(s_j+1/2/ρ_j+1/2,l){ψ_j+1,l - ψ_j,l} - (s_j-1/2/ρ_j-1/2,l){ψ_j,l - ψ_j-1,l}] + (1/s_j²)⋅(1/ΔΦ²)⋅[(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 = k₁ψ_j,l+1 + k₂ψ_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 = k₁ψ_j,l+1 + k₂ψ_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 R₁ on an initial guess, solve for Ψ₁,l in terms of Ψ_2,l and then find R₂ 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?

 

[phyzguy]

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.

 

[Daniel Sellers]

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.

 

[Fred Wright]

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.

 

[Daniel Sellers]

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!

 

[Daniel Sellers]

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 = ρ²c1 per the first substitution?

For the second set of substitutions, if we replace Ψ in the eq by Ψ' I can't get a₀ 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.