I am trying to solve a system of the form
where is a decreasing function of x. It's similar to a power law but not quite, so I expect to also be nearly a power law. I have a boundary condition , and . The solution I want will be the initial condition of a differential equation in time, so obviously I want to get this as right as possible.
My problem is this: the values of are defined on a grid with varying grid spacing and I am using a centered, three-point finite difference formula to represent derivatives at each point (except at the outer boundary, where I use a one-sided one instead). If I merely set up my array of y values the naive way (ie. I find that the result does not satisfy the differential equation using the finite difference derivatives. They start off in agreement, but by the time I get to the last few grid points they differ by four orders of magnitude.
Instead, I have tried solving it using the system of simultaneous equations defined by the finite difference formula:
where the come from the finite difference formula (setting :
This system is a tridiagonal matrix and is easily solved, I just have to be careful about the last row because it's a one-sided FD rather than a centered one but that's not much of a problem. Unfortunately, the result I get starts off well but finishes up zigzagging between two very different functions, neither of which resemble the approximate power law I expect. I've been tearing my hear out trying to figure out what's going on. I should be able to get an approximate solution to a well-behaved differential equation using long double precision in C, and it annoys me that I can't. Has anyone seen this kind of behaviour before? Is it typical of the kind of FD formula I am using, and if so how can I deal with it?
Alternatively, if you think my whole approach is stupid and you can suggest a better one, I am all ears.