I am trying to understand/develop a numerical code for solving PDEs. I started with a simple first order PDE, $$ \frac{\partial F(r,\theta)}{\partial r}+\frac{\partial F(r,\theta)}{\partial \theta} = 0 $$ Under the boundary condition that at maximum r $F(r_{max},\theta)=55$ (in my code $r_{max}=3$). Solving the PDF using backward finite difference method I came up with the following solution, $$ F(r_{i-1},\theta_i) = \frac{(h_r+h_\theta) F(r_{i},\theta_i) - h_r F(r_{i},\theta_{i-1})}{h_\theta} $$ Here, $h_r$ and $h_\theta$ are the step size for $r$ and $\theta$ respectively. I wrote the following Matlab code for for this solution,
hr = 0.1;
r = 0:hr:3;
hT = 0.1;
T = 0:hT:2*pi;
F = zeros(length(T), length(r));
F(:,length(r))=10;
for i=length(r):-1:2
for j=length(T):-1:2
F(j,i-1) = (((hr+hT)*F(j,i)) - (hr*F(j-1,i)))/hT;
end
end
surf(F)
Based on this I am getting the following surface as solution,
So here are my questions: 1) My solution based on the finite difference shows that having a spherically symmetric boundary condition (no $\theta$ dependence) means that the whole solution is spherically symmetric(?). Is it reflected in the surface plot? I am not sure. May be it is the cartesian coordinate frame that is restricting my imagination to see that this solution actually is spherically symmetric.
2) The solution has weird behavior. For example, if I change the step size in $r$ to a different value, like $h_r = 0.01$ then I get the following surface.
Notice the the scale of $F(r,\theta)$ has shifted by orders of magnitude. This is very odd in my opinion.
3) If I change the resolution of $\theta$ I get the following solution curve,
. Which is odd still because for some reason the structure has shifted and also the scale has changed.
So I am not able to reconcile all these behaviors and explain if this is expected. Unless the solution is highly non-linear so when I change the step size the scale changes (because more non-linear points are included which have extremely high value of the function) but that should be the case for such a simple differential equation, unless of course the PDE admits a solution which actually is this non-linear for the given boundary condition.
Adding a few more cases for reference,
$h_r = 0.1$ and $h_\theta=0.5$
$h_r = 0.5$ and $h_\theta=0.1$
I do not understand how is it possible that by changing the resolution of $r$ has effect on the theta resolution. Or may be
I am looking at it all wrong.
Look at your equation, and imagine you solve it in exact arithmetic. For $r=r_{max}$, it is constant in $\theta$. If it is a constant $C$ for $r=r_i$, $$ F(r_{i-1},\theta_i)= \frac{(h_r+h_\theta) F(r_{i},\theta_i) - h_r F(r_{i},\theta_{i-1})}{h_\theta} = \frac{(h_r+h_\theta) C - h_r C}{h_\theta}=C, $$ so all of your solution would be $C.$ The differences occur because you introduce rounding errors in every step and divide them by a small number, $h_\theta$. I'd imagine the following form is mathematically equivalent, but more stable: $$ F(r_{i-1},\theta_i)= F(r_{i},\theta_i)+h_r\,\frac{F(r_{i},\theta_i) - F(r_{i},\theta_{i-1})}{h_\theta} $$ I wonder what would be the physical meaning of such an equation in polar coordinates, though. It's easy to see that its general solution is $F(r,\theta)=f(r-\theta)$, meaning it's constant along the lines $r=\theta+c.$ Those lines are Archimedean spirals, and it's clear you get a singularity, an eddy, at $r=0.$