Successive substitution function in C

Question

I have to create a function that does successive substitution.

For example, I have a function F(X,Y) = (sqrt ((X-Y)/(X-1))) / tanh(sqrt((X-Y)/(X-1))) .

(X is a constant)

I have to enter and check values of Y (smaller than X) like this:

    Value 1: SET Y=0
         E1 = F(X,Y)
  then   E2 = F(X,E1)
  THEN   E3 = F(X,E2)
          .
          .
          .
          .
        E(n) = F(X,E(n-1))

and stop when (fabs(E(n)-E(n-1))<= 0.001); 

I tried to use a do-while loop but I didn't manage to make it work.

#include <math.h>
#include <stdio.h>
#include <stdlib.h>

#define F(Ei,E)   (sqrt (M*(Ei-E)/(Ei-1)))  /  tanh(sqrt(M*(Ei-E)/(Ei-1)))

#define Ei        1 + ((Db*CBoo)/(b*Da*CAi))
#define M         (Da*k2*CBoo)/(kl*kl)

int main(){

    double k2,Da,Db,b,h,CAi,CBoo,kl,E,i;
    double E1 ,E2 ; 
    kl =1;   
    k2 = 3;
    Da = 0.4;
    Db = 0.3;
    CAi = 3;
    CBoo = 5; 
    b = 2;
    E1=0;
    E2=0;
    E=0;
    i=0;

    do
    {
        E1 = E2;
        E2 = F(Ei,E1);

        printf("dd");
    }while(fabs(E2-E1)<= 0.0001);

    printf("\nanswer is %lf",E1);

    system("PAUSE");
}

The loop only does 1 round and the answer is always 0

Any ideas how to do this??

Thank you!

Answer 1

First, you need to change while(fabs(E2-E1) <= 0.001) to while(fabs(E2-E1) > 0.001) as was pointed out by MicroAlex.

The reason it still doesn't work is that the term sqrt(X-Y) restricts the domain of your function F(X,Y) to the interval [0, X] unless you are prepared to deal with complex numbers (which you are not, obviously). With your starting point Y0=0, you immediately get a new Y1=F(X,Y0) outside the domain of your function. This is why it breaks after one iteration. As a matter of fact you're lucky it breaks since the behaviour of the loop condition is essentially undefined.

So the real problem lies in your numerical approach. The problem doesn't seem to yield to fixed-point iteration. What are you trying to do, actually?

edit

The reason why fixed-point iteration doesn't work is probably that the function is not Lipschitz-continuous (due to the sqrt term, but I haven't checked). You might be better off with trying a different root solver, like bisection with a start interval [0, X]

Answer 2

You should use while(fabs(E2-E1) > 0.001) instead of while(fabs(E2-E1) <= 0.001)