Método de Muller en lenguaje C

#include
#include
#include "complejo.h"
#define CERO 0.00000001


fcomplex f(fcomplex x){
return Cadd(Cmul(x,Cmul(x,x)),Complex(-27.,0.));/*x^3 -27*/
}

fcomplex df(fcomplex x){
return RCmul(3.,Cmul(x,x));/*3x^2*/
}
/****************************Newton Rapson***********************************/

fcomplex Newton_Rapson(fcomplex p, float tol, int n){
int i=0;
fcomplex np;
np=p;
p = Csub(p, Cdiv(f(p),df(p)) );

while ((Cabs(Csub(np,p)) > tol) && (Cabs(f(p)) > CERO) && (i i++;
np=p;

if(Cabs(df(p))>CERO)
p =Csub(p,Cdiv(f(p),df(p)));

}
if (i==n)
printf("NEWTON HA EXEDIDO EL NUMERO DE ITERACIONES");
return(p);
}

/**********************calculo de la segunda derivada***********************/


fcomplex d2 (fcomplex aprox1, fcomplex aprox2,fcomplex aprox3)



{
fcomplex a,b,c;
if ( (Cabs(Csub(aprox2,aprox1)) > CERO) && (Cabs(Csub(aprox3,aprox2)) > CERO) ) {
a = Cdiv(Csub(f(aprox2),f(aprox1)),Csub(aprox2,aprox1));
b = Cdiv(Csub(f(aprox3),f(aprox2)),Csub(aprox3,aprox2));
c = RCmul(2.0,(Cdiv(Csub(a,b),Csub(aprox1,aprox3)))); }
else
printf("\nDIVISION POR CERO EN DERIVADA SEGUNDA\n");
return c;
}

/*definicion de la primera derivada discretizadas*/



fcomplex d1 (fcomplex aprox1, fcomplex aprox2, fcomplex aprox3)



{
fcomplex a,b,c;
if (Cabs(Csub(aprox3,aprox2)) > CERO )
a = Cdiv(Csub(f(aprox3),f(aprox2)),Csub(aprox3,aprox2));
else
printf("\nDIVISION POR CERO EN DERIVADA PRIMERA\n");
b = Cmul(RCmul(0.5,d2(aprox1,aprox2,aprox3)),Csub(aprox3,aprox2));
c = Cadd(a,b);
return c;
}




/***********************calculo de la raiz cuadrada*************************/

fcomplex Msqrt (fcomplex x1, fcomplex x2, fcomplex x3){

return Csqrt( Csub ( Cmul(d1(x1,x2,x3),d1(x1,x2,x3)) , Cmul( RCmul (2,d2(x1,x2,x3)), f(x3) ) ));
}


/***********************funcion Muller****************************/

fcomplex muller (fcomplex x0, fcomplex x1, fcomplex x2, float tol, int iter){

fcomplex raiz_pos, raiz_neg;
int i=0;

while ( ( Cabs(Csub(x1,x2)) > tol) && (iCERO) ){

if (Cabs(d2(x0,x1,x2)) > CERO){
raiz_pos = Cdiv( Cadd(RCmul(-1.,d1(x0,x1,x2)) ,Msqrt(x0,x1,x2)) ,d2(x0,x1,x2)) ;
raiz_neg = Cdiv( Csub(RCmul(-1.,d1(x0,x1,x2)) ,Msqrt(x0,x1,x2)) ,d2(x0,x1,x2)) ;

if (Cabs(raiz_pos) > Cabs(raiz_neg)){
x0=x1;
x1=x2;
x2=Cadd(raiz_neg,x2); }
else{
x0=x1;
x1=x2;
x2=Cadd(raiz_pos,x2); }
}
else{
printf("\nSOLUCION NEWTON RAPSON\n");
return Newton_Rapson(x0, tol, iter);}
i++;
}
if (i>=iter)
printf("\nMULLER HA EXEDIDO EL NUMERO DE ITERACIONES\n");

return x2;
}


/**************************************************************************************/

int main(){

fcomplex s;

s= muller ( Complex(0.01,0.), Complex(0.1,0.), Complex(1.,0.) ,0.000001, 1000);
printf("\nLa raiz buscada es: (%f,%f)\n\n", s.r,s.i);
return 0;
}

// complejo.h


#include
#include

#ifndef _FCOMPLEX_DECLARE_T_
typedef struct {float r,i;} fcomplex; /* type for a complex no */
#define _FCOMPLEX_DECLARE_T_
#endif /* _FCOMPLEX_DECLARE_T_ */

/* function prototypes */
/***********************/
fcomplex Cadd( fcomplex a, fcomplex b);
fcomplex Csub( fcomplex a, fcomplex b);
fcomplex Cmul( fcomplex a, fcomplex b);
fcomplex Complex( float re, float im);
fcomplex Conjg( fcomplex z);
fcomplex Cdiv( fcomplex a, fcomplex b);
float Cabs( fcomplex z);
fcomplex Csqrt( fcomplex z);
fcomplex RCmul( float x, fcomplex a);
fcomplex Cinv( fcomplex z);


// complejo.c

#include "complejo.h"

fcomplex Cadd(fcomplex a, fcomplex b)
{
fcomplex c;

c.r=a.r+b.r;
c.i=a.i+b.i;
return c;
}


fcomplex Csub(fcomplex a, fcomplex b)
{
fcomplex c;
c.r=a.r-b.r;
c.i=a.i-b.i;
return c;
}


fcomplex Cmul(fcomplex a, fcomplex b) /* Using only 3 multiplications! */
{
fcomplex c;
float t1,t2;

t1=a.r*b.r;
t2=a.i*b.i;

c.r=t1-t2;
c.i=(a.r+a.i)*(b.r+b.i)-(t1+t2);
return c;
}


fcomplex Complex(float re, float im)
{
fcomplex c;
c.r=re;
c.i=im;
return c;
}


fcomplex Conjg(fcomplex z)
{ fcomplex c;
c.r=z.r;
c.i = -z.i;
return c;
}

fcomplex Cdiv(fcomplex a, fcomplex b)
{
fcomplex c;
float r; float den;

if (fabs(b.r) >= fabs(b.i)) {
r=b.i/b.r;
den=b.r+r*b.i;
c.r=(a.r+r*a.i)/den;
c.i=(a.i-r*a.r)/den;
}
else {
r=b.r/b.i;
den=b.i+r*b.r;
c.r=(a.r*r+a.i)/den;
c.i=(a.i*r-a.r)/den;
}
return c;
}

float Cabs(fcomplex z)
{
float x,y,ans,temp;
x=fabs(z.r);
y=fabs(z.i);
if (x == 0.0)
ans=y;
else if (y == 0.0)
ans=x;
else if (x > y) {
temp=y/x;
ans=x*sqrt(1.0+temp*temp);
}
else {
temp=x/y;
ans=y*sqrt(1.0+temp*temp); }
return ans;
}


fcomplex Csqrt(fcomplex z)
{
fcomplex c;
float w;
if ((z.r == 0.0) && (z.i == 0.0)) {
c.r=0.0;
c.i=0.0;
}
else {
w = sqrt((sqrt(z.r*z.r + z.i*z.i) + fabs(z.r)) * 0.5);
if (z.r >= 0.0){
c.r=w;
c.i=z.i/(2.0*w);
}
else {
c.i=(z.i >= 0) ? w : -w;
c.r=z.i/(2.0*c.i);
}
}
return c;
}


fcomplex RCmul(float x, fcomplex a)
{
fcomplex c;
c.r=x*a.r;
c.i=x*a.i;
return c;
}


fcomplex Cinv( fcomplex z)
{
fcomplex c;
float s = 1.0 / (z.r*z.r + z.i*z.i);

c.r = z.r * s;
c.i = -z.i * s;
return c;
}


/* int main()
{
fcomplex a, b, c;

a = Complex( 10, 6 );
b = Complex( 4, 9 );

c = Cadd( a, b );
printf( "Suma: (%f,%f)\n", c.r, c.i );
c = Csub( a, b );
printf( "Resta: (%f,%f)\n", c.r, c.i );
c = Cmul( a, b );
printf( "Producto: (%f,%f)\n", c.r, c.i );
c = Cdiv( a, b );
printf( "Divisin: (%f,%f)\n", c.r, c.i );
} */