The Finite Difference method is based on discretizing the domain where our function is defined and then solve the equations by approximating the derivatives as differences of the value of our function at each point. THus far so good.
However, my problem is the following.
Let's say our domain is a unit square to simplify things, i.e $[0,1]^2.$
Imagine you are also given a function $f(x,y)$ such that
$$ f(x,y)=b, (x-0.5)^2+y(-0.5)^2< 0.25 \\ f(x,y)=a, otherwise $$
in essence, f is continuous and constant at all points of its domain except at the circumference. So it is basically a Heaviside stepfunction shifted upwards, but I'll use the above notation for clarity.
Now, in that domain, you want to solve the following equation: $$ \nabla \cdot (f.\vec{\nabla}u)=0 $$
where u(x,y) $I\!R^2 -> I\!R$ is the function we are interested in studying and solving for using the finite difference method.
My question is rather simple: how do we deal with the fact that $f$ is discontinuous at those points? Because if not for that, this would be really simple; but I confess I am stumped. I first thought of approximating f by a continuous function such as a fermi-dirac or even a polynomial function, but I feel there must be a simpler way.
Thank you for your time!
I think there is a simpler solution: remember the physics behind your problem. Your equation is derived from the absence of sources, for every closed domain, there flows as much in as out. If you consider a small rectangle where everything is twice differentiable, you'll get the equation as a limit. Along the line of discontinuity, you still have the condition that there flows as much into the outside as out of the inside, so $f\,\frac{\partial u}{\partial n}$ (derivative in normal direction) must be equal at both sides of that line. So you approximate your equation with second order differences in inner points of both regions, and with first order differences along the boundary, with $f=b$ inside your circle, and $f=a$ outside. I think you can easily fill in the details.