Results 1 to 3 of 3

Thread: Difficulties with finite difference formula

  1. #1
    Feb 2011

    Difficulties with finite difference formula

    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.

    Last edited by TheBigH; Aug 3rd 2011 at 12:17 PM. Reason: fixing typo
    Follow Math Help Forum on Facebook and Google+

  2. #2
    MHF Contributor
    Jester's Avatar
    Dec 2008
    Conway AR

    Re: Difficulties with finite difference formula

    The solution of your problem is

    $\displaystyle y = K - \int_0^x \dfrac{r(t)}{t^2}dt$

    Why not try using some the the usual techniques to numerically integrate the RHS (although I see a potential problem at $\displaystyle x = 0$)?

    BTW - what does $\displaystyle r(x)$ look like?
    Follow Math Help Forum on Facebook and Google+

  3. #3
    Feb 2011

    Re: Difficulties with finite difference formula

    I made a typo in my original post, sorry about that; it should read y(x_0) rather than y(0) so that the solution to the problem is
    $\displaystyle y=\int_x_0^x\frac{r(t)}{t^2}dt$, and the problem with zero is avoided.

    $\displaystyle r(x) \approx Cx^{-4/3}$, with the deviations from the power law less than about one part in a thousand.

    I have already tried some of the usual techniques for numerically integrating the RHS. This of course gets me a function that looks like what I expect, but when I apply my finite difference scheme to it I get errors that are small compared to the function itself but still a few orders of magnitude higher than I am able to tolerate. That's why I am doing this thing with the tridiagonal matrix- I figure that I should set up the y(x) array so that applying the FD formula to it gives me exactly -r(x)/x^2.

    I find that this works stunningly well in terms of accuracy. My errors are reduced to about a billionth of what they were, and well within what I can afford to put up with. The trouble is that the resulting function does not look like what I expect, particularly for high x values. There seem to be two distinct solutions that y(x) zigzags between, and I still get humungous errors at the inner and outer boundaries.

    By playing around with my code, I have found that I can suppress the worst of it by increasing the number of grid points and making the grid spacing very small near the lower boundary. But I'm still confused as to why it happens at all.
    Follow Math Help Forum on Facebook and Google+

Similar Math Help Forum Discussions

  1. Eigenvalues a finite difference formulas
    Posted in the Advanced Algebra Forum
    Replies: 0
    Last Post: Oct 6th 2011, 02:54 AM
  2. Replies: 2
    Last Post: Aug 12th 2011, 01:39 AM
  3. Finite Difference Problem
    Posted in the Advanced Applied Math Forum
    Replies: 1
    Last Post: May 21st 2009, 03:02 PM
  4. Explicit Finite difference
    Posted in the Differential Equations Forum
    Replies: 0
    Last Post: Apr 1st 2009, 07:12 AM
  5. difference between finite element/finite difference method
    Posted in the Advanced Applied Math Forum
    Replies: 0
    Last Post: Oct 2nd 2008, 11:03 AM

Search Tags

/mathhelpforum @mathhelpforum