Hi everyone,

I am trying to solve a system of the form

$\displaystyle \frac{dy}{dx}=-\frac{r(x)}{x^2}$

where $\displaystyle r(x)$ is a decreasing function of x. It's similar to a power law but not quite, so I expect $\displaystyle y(x)$ to also be nearly a power law. I have a boundary condition $\displaystyle y(x_0)=K$, and $\displaystyle x_0>0$. 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 $\displaystyle y(x)$ 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. $\displaystyle y(x_i) = y(x_{i-1}) - \frac{r(x_i)}{x^2_i}(x_i-x_{i-1}))$ 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:

$\displaystyle \left[\begin{array}{cccccc}1&0&0&0&0&\dots\\a_1&b_1&c_1& 0&0&\dots\\0&a_2&b_2&c_2&0&\dots\\0&0&a_3&b_3&c_3& \dots\\\vdots&\vdots&\vdots&\vdots&\vdots&\ddots\\ \end{array}\right]\left[\begin{array}{c}y_0\\y_1\\y_2\\y_3\\ \vdots\\\end{array}\right]=\left[\begin{array}{c}K\\r_1/x_1^2\\r_2/x_2^2\\r_3/x_3^2\\\vdots\\\end{array}\right]$

where the $\displaystyle a_i,b_i,c_i$ come from the finite difference formula (setting $\displaystyle h=x_i-x_{i-1},k=x_{i+1}-x_i$:

$\displaystyle \frac{dy}{dx}\approx \frac{-ky_{i-1}}{h(h+k)}+\frac{(k-h)y_i}{hk}+\frac{hy_{i+1}}{k(h+k)}$

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.

Thanks,

H.