Suma de dos vectores en lenguaje ADA

La suma de dos vectores (V1, V2) de igual tamaño, es un nuevo vector, V3, del mismo tamaño que V1 y V2, tal que, cada elemento de V3 es la suma de los que ocupan las posiciones correspondientes (primera, segunda, tercera,...; independientemente de los índices que designen dichas posiciones) de V1 y V2.

Implemente una función de tipo Vector_1, llamada Suma_Vectores, con dos parámetros de tipo Vector_1. El resultado será la suma de los dos parámetros de la función. Dado que dos vectores sólo se pueden sumar cuando tienen el mismo número de elementos, la función lanzará una excepción Constraint_Error con el mensaje "Vectores incompatibles" cuando el tamaño de los parámetros de entrada sea distinto. El rango del índice del resultado deberá empezar siempre en 1, independientemente de cuáles fueran los rangos de índices de los parámetros de la función.

with Arrays; use Arrays;
function Suma_Vectores (V1, V2 : Vector_1) return Vector_1 is
suma : Vector_1 (1 .. V1'Length);
begin
if V1'Length = V2'Length then
for I in 1 .. V1'Length loop
suma (I) := V1 (V1'First + I - 1) + V2 (V2'First + I - 1);
end loop;
else
raise Constraint_Error with "Vectores incompatibles";
end if;
return suma;
end Suma_Vectores;

package Arrays is
type Vector_1 is array (Integer range <>) of Integer;
type Vector_2 is array (Positive range <>) of Integer;
type Matriz_1 is array (Integer range <>, Integer range <>) of Integer;
end Arrays;

Fila de mayor peso en una matriz en lenguaje ADA

Se define el peso de una fila de una matriz de enteros como la suma de todos sus elementos.

Implemente una función de tipo Integer, llamada Fila_Mayor_Peso, con un parámetro de tipo Matriz_1. La función devolverá el índice de la fila de mayor peso de la matriz pasada como parámetro. Si hay varias filas con el mismo peso se devolverá el índice menor.

with Arrays; use Arrays;
function Fila_Mayor_Peso (Matriz : Matriz_1) return Integer is
Fila_Maxima : Integer;
Indice_Fila_Maxima : Integer;
aux : Integer := Integer'First;
begin
for I in Matriz'Range (1) loop
Fila_Maxima :=1;
for J in Matriz'Range (2) loop
Fila_Maxima := Fila_Maxima + Matriz (I, J);
end loop;

if Fila_Maxima > aux then
Indice_Fila_Maxima := I;
aux := Fila_Maxima;
end if;
end loop;
return Indice_Fila_Maxima;
end Fila_Mayor_Peso;

package Arrays is
type Vector_1 is array (Integer range <>) of Integer;
type Vector_2 is array (Positive range <>) of Integer;
type Matriz_1 is array (Integer range <>, Integer range <>) of Integer;
end Arrays;

Apareciones de un elemento en un vector en lenguaje ADA

Implemente una función de tipo Natural, llamada Apariciones, con dos parámetros: el primero de tipo Vector_1 y el segundo de tipo Integer. La función debe devolver el número de veces que el valor representado por el segundo parámetro aparece en el array representado por el primero.


with Arrays; use Arrays;
function Apariciones (Vector : Vector_1; Valor : Integer) return Natural is
contador : Natural := 0;
begin
for i in Vector'Range loop
if Vector (i) = Valor then
contador := contador + 1;
end if;
end loop;
return contador;
end Apariciones;

package Arrays is
type Vector_1 is array (Integer range <>) of Integer;
type Vector_2 is array (Positive range <>) of Integer;
type Matriz_1 is array (Integer range <>, Integer range <>) of Integer;
end Arrays;

Número de positivos en un vector en lenguaje ADA

Implemente una función de tipo Natural, llamada Positivos, con un parámetro de tipo Vector_1, que devuelva el número de valores positivos contenidos en el array representado por dicho parámetro.

with Arrays; use Arrays;
function positivos (Vector : Vector_1) return Natural is
contador : Natural := 0;
begin
for I in Vector'Range loop
if Vector (I) > 0 then
contador := contador + 1;
end if;
end loop;
return contador;
end positivos;

package Arrays is
type Vector_1 is array (Integer range <>) of Integer;
type Vector_2 is array (Positive range <>) of Integer;
type Matriz_1 is array (Integer range <>, Integer range <>) of Integer;
end Arrays;

Practica en Lenguaje ADA: Implementa el préstamo Francés y Americano. También una cola FIFO, para las existencias en almacén.

with unchecked_deallocation,Text_Io,ada.integer_text_io,Ada.Float_Text_IO;
use Text_Io,ada.integer_text_io,Ada.Float_Text_IO;


procedure trabadministracion is
type nodo_lista;
type pnodo_lista is access nodo_lista;
type nodo_lista is record --lista doblemente encadenada
cantidad:integer;
precio:float;
siguiente:pnodo_lista;
anterior:pnodo_lista;
end record;
principio,final:pnodo_lista;--punteros al nodo ultimo y primero
---variables de los prestamos--
cantidad,cant,opcion:integer;
aux,precio:float;
---variables del lifo--
n,cont,s:integer;
g,anual,interes,i,c,f,a:float;

procedure Liberar_pnodo is new Unchecked_Deallocation(nodo_lista,Pnodo_lista);
--libera nodos del tipo pnodo_lista


procedure insertarfinal(cantidad: in integer;Precio: in float;principio,final: in out pnodo_lista) is
--inserta un nuevo nodo al final de la lista
recorrer:pnodo_lista;
begin
recorrer:=final;
if principio=null then --en caso de lista vacia
principio:=new nodo_lista;
principio.cantidad:=cantidad;
principio.precio:=precio;
final:=principio;
else --si la lista no esta vacia
final:=new nodo_lista;
final.cantidad:=cantidad;
final.precio:=precio;
recorrer.siguiente:=final;
final.anterior:=recorrer;
end if;
end insertarfinal;

procedure extraer(cant: in out integer;final: in out pnodo_lista)is
aux,aux2:pnodo_lista;
Cant_total:integer:=0 ;
begin
aux2:=final;
while aux2/=null loop
cant_total:=cant_total+aux2.cantidad;
aux2:=aux2.anterior;
end loop;
if cant>cant_total then
new_line;
put_line("La Cantidad Es Superior A La Actual");
else
if cant=final.cantidad then
aux:=final.anterior;
if principio/=final then
aux.siguiente:=null;
else
principio:=aux;
end if;
liberar_pnodo(final);
final:=aux;

else
if cant>final.cantidad then
cant:=cant-final.cantidad;
aux:=final.anterior;
liberar_pnodo(final);
aux.siguiente:=null;
final:=aux;
extraer(cant,final);
else
final.cantidad:=final.cantidad-cant;
end if;
end if;
end if;
end extraer;

procedure verpantalla(principio,final:pnodo_lista) is
--visualiza en pantalla la lista
recorrer:pnodo_lista;
begin
recorrer:=principio;
if not(recorrer=null) then --si la lista no esta vacia
--recorre toda la lista y pone el campo info en pantalla
while not(recorrer=null) loop
put("cantidad:");put(recorrer.cantidad);
put("precio:");put(recorrer.precio,fore=>5,aft=>2,exp=>0); --te pasa de la notacion cientifica a float
new_line;
recorrer:=recorrer.siguiente;
end loop;


else --si la lista esta vacia llama la excepcion
new_line;
put_line("No Hay Existencias");
end if;
end verpantalla;

function anualidad(c,i:float;n:integer)return float is
f,a,i2:float;
begin
i2:=i/100.0;
f:=1.0+i2;
f:=f**(-n);
f:=1.0-f;
f:=f/i2;
a:=c/f;
return (a);
end anualidad;

function cinteres(c,i,a:float;n,s:integer)return float is
f,i2,cinteres:float;
s2,n2:integer;
begin
s2:=s-1;
n2:=n-s2;
i2:=i/100.0;
f:=1.0+i2;
f:=f**(-n2);
f:=1.0-f;
f:=f/i2;
cinteres:=f*(anualidad(c,i,n))*i2;

return cinteres ;
end cinteres;

function camortizacion(c,i,a:float;n,s:integer)return float is
f,i2,i3,cinteres,a2,a3:float;
s2,n2:integer;
begin
i2:=i/100.0;
i3:=c*i2;
a2:=(anualidad(c,i,n))-i3;
s2:=s-1;
a3:=(1.0+i2)**s2;
a3:=a2*a3;

return a3 ;
end camortizacion;

function saldo(c,i,a:float;n,s:integer)return float is
f,i2,a2:float;
s2:integer;
begin
s2:=n-s;
i2:=i/100.0;
f:=1.0+i2;
f:=f**(-s2);
f:=1.0-f;
f:=f/i2;
a2:=(anualidad(c,i,n))*f;
return a2 ;
end saldo;
begin
loop
--begin
new_line;
put_line("*********************** Bienvenido A Este Programa *************************");
put_line("************** Administracion De Empresas II **************");
put_line("******* *******");
put_line("******* 1.= Terminar EL programa *******");
put_line("******* 2.= Realizado por *******");
put_line("******* 3.= Lifo *******");
put_line("******* 4.= Prestamos *******");
put_line("******* *******");
put_line("****************************************************************************");
new_line;
Put("Teclee su opcion");
get(opcion);
new_line;
case Opcion is
when 1=>
put_line("gracias por usar este programa");
exit;
when 2=>
new_line;
Put_line(" ***** Diseñadores Del Econominator ******");
new_line;
Put_line(" ** Jose Luis Cabrera Marrero ** ");
new_line;
Put_line(" ** Airam Caballero Perez **");
new_line;
Put_line(" ** Eduardo Hidalgo Luque **");
new_line;
Put_line(" ** Sergio Romero Sanchez **");
New_Line;
When 3=>
loop
new_line;
put_line("******************************* LIFO ********************************");
put_line("************** Administracion De Empresas II *************");
put_line("******* *******");
put_line("******* 1.= menu principal *******");
put_line("******* 2.= entradas *******");
put_line("******* 3.= salidas *******");
put_line("******* 4.= existencias actuales *******");
put_line("******* *******");
put_line("****************************************************************************");
new_line;
Put("Teclee su opcion");
get(opcion);
new_line;
case opcion is
when 1=>
EXIT;
when 2=>
put_line("-ENTRADAS-");
put("indique la cantidad de mercancia entrante:");
get(cantidad);
new_line;
put("precio al que entra la mercancia:");
get(precio);
insertarfinal(cantidad,precio,principio,final);
new_line;
put_line("¡¡ la entrada ha sido realizada ¡¡");
when 3=>
put_line("cantidad que va a salir del almacen");
get(cant);
extraer(cant,final);
new_line;

When 4=>
put_line("las existencia actuales son:");
verpantalla(principio,final);
When Others=>
Put_line("La Opcion Elegida Es Incorrecta");
New_Line;
New_Line;

End case;
end loop; ----fin del menu lifo
new_line;
when 4=>
loop
new_line;
put_line("******************************* PRESTAMOS ***************************");
put_line("************** Administracion De Empresas II *************");
put_line("******* *******");
put_line("******* 1.= menu principal *******");
put_line("******* 2.= Sistema Americano *******");
put_line("******* 3.= Sistema Frances *******");
put_line("******* *******");
put_line("****************************************************************************");
new_line;
Put("Teclee su opcion");
get(opcion);
new_line;
case opcion is
when 1=>
EXIT;
when 2=>
loop
new_line;
put_line("***************************** PRESTAMOS **************************");
put_line("**** --SISTEMA AMERICANO-- ****");
put_line("**** ****");
put_line("**** 1.= menu principal ****");
put_line("**** 2.= cuantia constante de interes a pagar en cada ano del prestamo ****");
put_line("**** 3.= cuantia constante que hay que imponer en el fondo ****");
put_line("**** ****");
put_line("*******************************************************************************");
new_line;
Put("Teclee su opcion");
get(opcion);
new_line;
case opcion is
when 1=>
exit;
when 2=>

put_line("introduzca el capital del prestamo");
get(c);
put_line("introduzca el interes del prestamo anual en %");
get(i);


begin
aux:=i/100.0;
interes:=c*aux;
new_line;
put_line("la cuantia de interes es :"); put(interes,fore=>5,aft=>3,exp=>0);

end;
when 3=>
put_line("introduzca el capital del prestamo");
get(c);
put_line("introduzca las imposiciones anules en %");
get(i);
put_line("introduzca la duracion del prestamo en años");
get(n);
begin
f:=i/100.0;
f:=(((1.0+f)**n)-1.0)/f;
cont:=0;
anual:=c/f;
new_line;
put_line("la cuantia es:"); put(anual,fore=>5,aft=>3,exp=>0);


end;
when others=>
put_line("La Opcion Es Incorrecta");
trabadministracion;
end case;
end loop; ---fin del menu sistema americano
when 3=>
put_line("introduzca el capital del prestamo");
get(c);
put_line("introduzca el interes del prestamo anual en %");
get(i);
put_line("introduzca la duracion del prestamo en años");
get(n);
loop

new_line;
put_line("**************************** PRESTAMOS ********************************");
put_line("************** ---SISTEMA FRANCES--- *************");
put_line("******* *******");
put_line("******* 1.= salir *******");
put_line("******* 2.= calcular la anualidad que amortiza el prestamo *******");
put_line("******* 3.= cuota de amortizacion de un periodo en concreto *******");
put_line("******* 4.= cuota de interes de un periodo concreto *******");
put_line("******* 5.= saldo del prestamo en un año concreto *******");
put_line("******* *******");
put_line("****************************************************************************");
new_line;
Put("Teclee su opcion");
get(opcion);
new_line;
case opcion is
when 1=>
exit;
when 2=>
put(anualidad(c,i,n),fore=>5,aft=>2,exp=>0);
when 3=>
put_line("introduzca el periodo que desea obtener:");
get(s);
if s>n then
put_line("el periodo seleccionado es mayor que los periodos de amortizacion del prestamo");
exit;
end if;
put_line("la cuota de amortizacion es:");
put(camortizacion(c,i,a,n,s),fore=>5,aft=>2,exp=>0);

when 4=>
put_line("introduzca el periodo que desea obtener:");
get(s);
if s>n then
put_line("el periodo seleccionado es mayor que los periodos de amortizacion del prestamo");
exit;
end if;
put_line("la cuota de interes es:");
put(cinteres(c,i,a,n,s),fore=>5,aft=>2,exp=>0);

when 5=>
put_line("introduzca el periodo que desea obtener:");
get(s);
if s>n then
put_line("el periodo seleccionado es mayor que los periodos de amortizacion del prestamo");
exit;
end if;
put_line("el saldo del prestamo es:");
put(saldo(c,i,a,n,s),fore=>5,aft=>2,exp=>0);
when others=>
put_line("La Opcion Es Incorrecta");
trabadministracion;


end case;
end loop; ----------fin del menu del prestamo fraces
when others=>
put_line("La Opcion Es Incorrecta");
trabadministracion;

end case;
end loop; ---fin del menu prestamos


when others=>
put_line("La Opcion Es Incorrecta");
trabadministracion;
end case;
--Creamos Las Excepciones Necesarias.").
Exception
When Data_error =>
new_line;
put_line("Los Parametros Introducidos Son Incorrectos");
skip_line;
new_line;
When others=>
Put_Line("Error Desconocido");
New_Line;
end;
end loop; ---fin menu principal
End trabadministracion;

Número menor en doble precisión en lenguaje C

#include
int main()
{
double min,num; /*Min; guardará el valor del mayor número real, y N el valor más grande de la mantisa*/
int e;
min=1.;
e=0;
num=1.;
while (min>min/2)
{
num=min;
min=min*0.5;
e=e-1;
}

printf("El menor número en simple presición tendrá como exponente %d y será el %e \n", e, num);
return 0;
}

Número menor en simple precisión en lenguaje C

#include
int main()
{
float min, min_aux; /* Min_aux; guardará el valor del menor número positivo */
int e,i;
min=1.;
e=0;
i=0;
while (min>min/2) /* hallamos el menor exponente y a la vez calculamos 2 elevado al menor exponente */
{
min_aux=min;
min=min/2;
e=e-1;
}
/*e=e+1;*/
min_aux*= 1./2; /* multiplicamos 2 elevado al menor exponente (min_aux) por la menor mantisa (1/2) */

printf("El menor número en simple presición tendrá como exponente %d y será el %e \n", e, min_aux);
return 0;
}

Número mayor en simple precisión en lenguaje C

#include
int main()
{
float max,n,res,aux;/*Max; guardará el valor del mayor número real, y N el valor más grande de la mantisa*/
int e,i;
max=1.;
e=0;
n=0.;
res=1.;
aux=1.;
while (max<2*max) /* aquí hallamos el exponente máximo,multiplicando 1 por 2 hasta que sea infinito*/
{

max=2*max;
e=e+1;
}
for (i=0; i<24; i++) /* hacemos la suma de la mayor mantisa (1/2 + 1/4...+1/2elev 24)*/
{
aux=aux*(1./2);
n=n+aux;
}
for (i=0; i {
res=res*2;
}
res=res*n*2;
printf("El mayor número en simple precisión tendrá como exponente %d y será %e \n", e,res );
return 0;
}

Número mayor en doble precisión en lenguaje C

#include
int main()
{
double max,n,res,aux; /*Max; guardará el valor del mayor número real, y N el valor más grande de la mantisa*/
int e,i;
max=1.;
e=0;
n=0.;
res=1.;
aux=1.;
while (max<2*max)
{
max=2*max;
e=e+1;
}
for (i=0; i<53; i++)
{
aux=aux*(1./2);
n=n+aux;
}
for (i=0; i {
res=res*2.;
}
res=res*n*2.;
printf("El mayor número en doble precisión tendrá como exponente %d y será %e \n", e, res );
return 0;
}

Método Newton Raphson en lenguaje C (solo)

#include
#include
#define CERO 0.000000001

/*Funcion que queremos aproximar*/
float f(float x){
return (x*x*x*x*x)+(x*x)-9;
}

/*Derivada de la función a aproximar*/
float df(float x){
return (5*(x*x*x*x)+2*x);
}

float Newton_Rapson(float p, float tol, int n){
int i=0;
float np;
np=p;

if(fabs(df(p))>CERO)
p = p-(f(p)/df(p));

while ((fabs(np-p) > tol) && (fabs(f(p)) > CERO) && (i i++;
np=p;
if(fabs(df(p))>CERO)
p =p-(f(p)/df(p));
}

if (i==n)
printf("\nSe alcanzo el numero maximo de iteraciones.\n");

return(p);
}

int main()
{
float sol;
sol=Newton_Rapson(1.5, 0.001, 50);
printf("\nLa solucion a Newton-Raphson es la raiz: %f\n", sol);
return 0;
system ("pause");
}

Método de Simpson en lenguaje C

#include
#include
#include


double potencia(double x, int grado) /*Calcula el cuadrado de cada uno de los componentes del vector*/
{
double resultado=1.;
int i;
if(grado==0)
return resultado;
else
{
for(i=0;i resultado = resultado*x;

}
return resultado;
}

double f(double *coef,double x,int grado)/*Aqui se calcula el valor de la funcion*/
{
int i;
double res = 0.;
for(i=0;i<=grado;i++)
res = res + coef[i]*potencia(x,i);
return res;
}


double simpson(double *coef,int grado,double x,double x1)
{
int i;
double n=100,lim,xm,simp=0.;
lim = x1-x; /*Limsup-Liminf*/
xm = lim/n;/*Lo dividimos para dividir el integral en subintegrales*/

for (i=0;i simp = simp + (f(coef,x+i*xm,grado)+f(coef,x+(i+1)*xm,grado)+4*f(coef,x+(i+0.5)*xm,grado))/6*xm;
return simp;
}
int main(){
double res, *vec;
double linf,lsup,com;
int i,tam=0;
printf("\t---------------------------------------------\n");
printf("\t| METODO DE SIMPSON |\n");
printf("\t---------------------------------------------\n");
printf("\nIntroduzca el grado del polinomio:");
scanf("%d",&tam);
vec=(double*)malloc(tam*sizeof(double));
for (i=tam;i>=0;i--)
{
printf("Introduzca el valor del coeficiente de x%d: ",i);
scanf("%lf",&com);
vec[i]=com;
}
printf("\nIntroduzca los limites del integral:");
printf("\nLimite inferior (xk):");
scanf("%lf",&linf);
printf("Limite superior (xk+1):");
scanf("%lf",&lsup);
res = simpson(vec,tam,linf,lsup);

printf("\n El Resultado de la integral es: %lf",res);
return 0;
}

Calculo de Secante en lenguaje C

#include
#include
#define CERO 0.000001

float f(float x){
return (x*x)-2.0;
}

float df(float x0,float x1){


return ((f(x1)-f(x0) )/ (x1-x0) );
}
/*****************************************************/
float Secante(float p0, float p1 , float tol, int n){
int i=0;
float p,q0,q1;

q0=f(p0);
q1=f(p1);

while ((fabs(p1-p0) > tol) && (fabs(f(p)) > CERO)) {
i++;

p = p1- q1*(p1-p0)/(q1-q0);

if( fabs(p-p1) > CERO )
return p;

p0=p1;
q0=q1;
p1=p;
q1=f(p);


}
printf("el n de iteraciones es %d\n", i);
return(p);
}

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

int main()
{
printf("La solucion es: %.16e\n\n", Secante(-2.0 , 2.0, 0.000001,10));
return 0;
}

Método Regula Falsi en lenguaje C

#include
#include

#define Cero 0.00001
float f(float x){
return ((x*x*x)-2);
}

float Regual_falsi (float a,float b, float tol){
float c = a - ( (b-a)/ (f(b)-f(a)) )*f(a);
int i=0;

while ( (f(c) != Cero) && ( fabs(b-a) > tol )&& i<10 ){


if (f(a) * f(c) < Cero)
b=c;
else
a=c;

c = a - ( (b-a)/ (f(b)-f(a)) )*f(a);
i++;}

return c;
}
int main(){

printf("el valor de la raiz es: %f \n", Regual_falsi (-2.0, 2.0, 0.000001 ) );
return 0;
}

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 );
} */

Método de Relajación en lenguaje C

#include
#include
#include

double** LeerDatos(char *nombrefichero,int *dimension)
{
double **matriz;
int i,j;
FILE *f;
if(f=fopen(nombrefichero,"r"),!f)
{
printf("Problema con la lectura del fichero\n");
return matriz;
}
fscanf(f,"%d\n",&(*dimension));
/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
{
matriz[i]=(double*)malloc((*dimension+1)*sizeof(double));
}
/* Lectura de elementos de la matriz*/
for(i=0;i<*dimension;i++)
{
for(j=0;j<*dimension;j++)
{
fscanf(f,"%lf\n",&(matriz[i][j]));
}
}
/* Leemos el vector de terminos independientes */
for(i=0;i<*dimension;i++)
{
fscanf(f,"%lf\n",&(matriz[i][*dimension]));
}
fclose(f);
return matriz;
}

/* Procedimiento para mostrar matrices */
void muestra(double **matriz,int Tam)
{
int i,j;
printf("\n\t");
for(i=0;i {
for(j=0;j {
printf("%lf ",matriz[i][j]);
}
printf("\n\t");
}
printf("\n");
}
/* Procedimiento que multiplica matrices */
double** mulmat(double **A, double **B,int dimension)
{
int i,j,k;
double **C=(double**)malloc(dimension*sizeof(double*));

for(i=0;i {
C[i]=(double*)malloc((dimension)*sizeof(double));
}
for(i=0;i {
for(j=0;j {
C[i][j]=0;
}
}
for (i=0;i {
for (j=0;j {
for (k=0;k {
C[i][j]=C[i][j]+(A[i][k]*B[k][j]);
}
}
}
return (C);
}
/* funcion que construye la matriz de Jacobi Mj*/
double** MJ(double **matriz,int dimension)

{
int i,j;
double **mj,**L_U,**Dinv;
/* Tomamos memoria */
L_U=(double**)malloc(dimension*sizeof(double*));
Dinv=(double**)malloc(dimension*sizeof(double*));
for(i=0;i {
L_U[i]=(double*)malloc((dimension)*sizeof(double));
Dinv[i]=(double*)malloc((dimension)*sizeof(double));
}
/* calculamos inv(D), y L_U=-L-U*/
for(i=0;i {
for(j=0;j {
if (j==i)
{
Dinv[i][j]=1/matriz[i][j];
L_U[i][j]=0;
}
else
{
Dinv[i][j]=0;
L_U[i][j]=(-1)*matriz[i][j];
}
}
}
/* mj=inv(d)*(-L-U)*/
mj=mulmat(Dinv,L_U,dimension);
free(L_U);
free(Dinv);
return (mj);
}
double norma(double *vec,int dimension,int Tnorma)
{
int i,cero=0.00000000001;
double temp=0;
if (Tnorma==1)
{/*Norma ecuclidea*/
for(i=0;i {
temp+=vec[i]*vec[i];
}
return sqrt(temp);
}
else
{/*Norma infinita*/
temp=fabs(vec[0]);
for (i=1;i {
if ((fabs(vec[i])-fabs(temp)) > cero)
temp=fabs(vec[i]);
}
return temp;
}
}

/* Devuelve el signo del autovalor*/
double sign(double *v1, double *v2, int dimension)
{
int i;
double temp=0;

for (i=0; i {
temp = temp + v1[i] * v2[i];
}
if (temp>=0)
return (1);
else
return (-1);
}

/* Multiplica una matriz por un vector*/
double* MatxVec(double** matriz,double* vec,int dimension)
{
int i,j;
double *sol=(double*)malloc(dimension*sizeof(double));

for (i=0;i {
sol[i]=0;
for (j=0;j {
sol[i]= sol[i] + matriz[i][j]*vec[j];
}
}
return sol;
}

double Potencia(double **matriz,int dimension,double *u_n,double tolerancia,int Nmax)
{
double* temp=(double*)malloc(dimension*sizeof(double));
double* u_n_1=(double*)malloc(dimension*sizeof(double));
double Norma,Norma_Anterior=0,solucion;
int n=1,j;

Norma=norma(u_n,dimension,1);
while((n<=Nmax)&&(fabs(Norma-Norma_Anterior)>tolerancia))
{
/*Guardamos la iteracion anterior*/
for (j=0;j {
u_n_1[j]=u_n[j];
}
/*Dividimos entre la norma la iteracion actual*/
for (j=0;j {
u_n[j]=u_n[j]/Norma;
}
/* Y resolvemos la iteracion siguiente*/
temp = MatxVec(matriz,u_n,dimension);
free(u_n);
u_n=temp;
Norma_Anterior=Norma;
/* Hallamos el signo del autovalor maximo*/
solucion=Norma*sign(u_n, u_n_1, dimension);
n=n+1;
Norma=norma(u_n,dimension,1);
}
if (n>Nmax)
printf("Se ha superado el numero de iteraciones maximo.\n");
free(u_n_1);
free(temp);
return (solucion);
}

void relajacion(double **matriz,double *b,double *u_n,int dimension,int Nmax,double w,double Tol)
{
double dif,temp,maxdif=1.;
double *anterior;
int i,j,k,n=1;

anterior=(double*)malloc(dimension*sizeof(double));
while((n<=Nmax)&&(maxdif>Tol)) /*Bucle para calcular las aproximaciones hasta el valor numero de iteraciones*/
{
printf(" u%d[ ",n);
for (i=0;i {
printf(" %lf",u_n[i]);
anterior[i]=u_n[i]; /*guardamos en el vector anterior la posicion anterior del vector de aproximacion*/
}
printf("]\n");

/*Calculamos la siguiente iteracion*/
for (i=0;i {
temp=0;
for (j=0;j {
if (i!=j)
{
temp=temp-(matriz[i][j]*u_n[j]);
}
}
u_n[i]=w*( (temp + b[i]) / matriz[i][i])+((1-w)*u_n[i]); /*Aqui se almacena la siguiente aprox obtenida*/
}
/*Calculamos la mayor diferencia para el calculo del error*/
maxdif=fabs(u_n[0]-anterior[0]);
for (k=1;k {
dif=fabs(u_n[k]-anterior[k]);
if (dif>maxdif)
{
maxdif=dif;
}
}
n++;
}
if (n>Nmax)
printf("Se ha superado el numero de iteraciones maximo.\n");
free(anterior);
}


int main()
{
int dimension,Nmax,i,j;
double Tol,wopt;
double **matriz,**Mj,*u,*b;

printf("\t---------------------------------------------\n");
printf("\t| METODO DE RELAJACION |\n");
printf("\t---------------------------------------------\n");
printf("\n Introduzca tolerancia: ");
scanf("%lf",&Tol);
printf(" Introduzca el numero maximo de iteraciones: ");
scanf("%i",&Nmax);

matriz=LeerDatos("matriz4.dat",&dimension);
u=(double*)malloc(dimension*sizeof(double));/*Vector de aproximaciones iniciales*/
b=(double*)malloc(dimension*sizeof(double));/*Vector de terminos independientes*/
printf(" Terminos independientes: ");
printf("\n b[");
for(i=0;i {
b[i]=matriz[i][dimension];
printf(" %lf",b[i]);
}
printf(" ]\n");
printf("\n La matriz de coeficientes es: \n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",matriz[i][j]);
}
printf("}\n");
}

/*Aplicamos el Metodo de la Potencia para despues calcular el Wopt.\n");*/
printf("\n Calculamos Mj=ind(D)(-L-U): \n");
Mj=MJ(matriz,dimension); /* Calculamos Mj=inv(D)(-L-U)*/
printf(" La matriz Mj da: \n");
muestra(Mj,dimension);
printf(" Introduce el vector de aproximacion inicial para el metodo de la potencia:\n");
for (i=0;i {
printf(" U[%d]=",i);
scanf("%lf",&(u[i]));
}

wopt=Potencia(Mj,dimension,u,Tol,Nmax); /* Calculamos el autovalor maximo*/

wopt=(2/(1+sqrt(1-(wopt*wopt))));/* Calculamos el w optimo*/

printf(" El Wopt es: %lf \n",wopt);
printf(" Aplicamos Relajacion.\n");
printf(" Introduce el vector de aproximacion inicial para el metodo de Relajacion:\n");

u=(double*)malloc(dimension*sizeof(double));
for (i=0;i {
printf(" U[%d]=",i);
scanf("%lf",&(u[i]));
}
/* Aplicamos relajacion */
relajacion(matriz,b,u,dimension,Nmax,wopt,Tol); /*procedure q nos devuelve de salida u*/
printf("\n El resultado[");
for(i=0;i {
printf(" %lf",u[i]);
}
printf(" ]\n");
return 0;
}

Método de Jacobi Mj_Cj en lenguaje C

#include
#include
#include

double** LeerDatos(char *nombrefichero,int *dimension)
{
double **matriz;
int i,j;
FILE *f;

if(f=fopen(nombrefichero,"r"),!f){
printf("Problema con la lectura del fichero\n");
return matriz;
}
fscanf(f,"%d\n",&(*dimension));


/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
matriz[i]=(double*)malloc((*dimension+1)*sizeof(double));


for(i=0;i<*dimension;i++)
for(j=0;j<*dimension;j++)
fscanf(f,"%lf\n",&(matriz[i][j]));


for(i=0;i<*dimension;i++)
fscanf(f,"%lf\n",&(matriz[i][*dimension]));

fclose(f);
return matriz;
}

void liberaMatriz(double **matriz,int dimension){
int i;

for(i=0;i free(matriz[i]);
free(matriz);
}
double **multiplica(double **A,double **B,int dimension){
int i,j,k;
double **C;
C=(double**)malloc(dimension*sizeof(double*));
for(i=0;i C[i]=(double*)malloc((dimension+1)*sizeof(double));

for (i=0;i for (j=0;j for (k=0;k C[i][j]=C[i][j]+A[i][k]*B[k][j];
return C;
}
double *matrizPorVector(double **matriz,double *vector,int dimension){
double *resultado,a,b,mult;
int i,j;

resultado=(double*)malloc(dimension*sizeof(double));

for(i=0;i mult=0;
for(j=0;j mult+=(matriz[i][j]*vector[j]);
resultado[i]=mult;
}

return resultado;
}

double *suma_vector(double *A,double *B,int dimension){
int i;
double *C;

C=(double*)malloc(dimension*sizeof(double));

for (i=0;i C[i] = A[i]+B[i];

return C;
}

void jacobi2 (double **m, double *v, double *aprox, int dimension, int iter){
double **D1,**LR, *C;
int i,j;


D1=(double**)malloc(dimension*sizeof(double*));
for(i=0;i D1[i]=(double*)malloc((dimension)*sizeof(double));

LR=(double**)malloc(dimension*sizeof(double*));
for(i=0;i LR[i]=(double*)malloc((dimension)*sizeof(double));


/*HALLAMOS LA MATRIZ D INVERSA*/

for (i=0;i for(j=0;j if (j==i)
D1[i][j] = 1./m[i][j];
else
D1[i][j] = 0.;



/*HALLAMOS LA MATRIZ -L-U*/

for (i=0;i for(j=0;j if (j!=i)
LR[i][j] = -m[i][j];
else
LR[i][j] = 0.;


/*OBTENEMOS LA MATRIZ DE JACOBI Mj*/
liberaMatriz(m,dimension);
m = multiplica(D1, LR, dimension);
printf("\n La matriz Mj da como resultado :\n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",m[i][j]);
}
printf("}\n");
}


/*CALCULAMOS EL VECTOR Cj*/
C=matrizPorVector(D1,v,dimension);

printf("\n El vector C da como resultado :\n\n");
for (j=0;j {
printf("%4.3f ",C[j]);
}

for(i=0;i aprox = matrizPorVector(m,aprox,dimension);
aprox = suma_vector(aprox, C, dimension);
}


printf("\n\nEl vector solucion es:\n");
for(i=0;i printf("\nvector Un[%d] = %lf",i,aprox[i]);


}

int main(){
double **matriz;
double *v, *aprox, aux;
int dimension, iter, i, j;

matriz=LeerDatos("matriz10.dat",&dimension);

printf("\t---------------------------------------------\n");
printf("\t| METODO DE GAUSS-JACOBI_MC |\n");
printf("\t---------------------------------------------\n");

v = (double*)malloc(dimension*sizeof(double));
aprox = (double*)malloc(dimension*sizeof(double));

printf("\n La matriz de coeficientes es: \n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",matriz[i][j]);
}
printf("}\n");
}

for (i=0;i v[i] = matriz[i][dimension];

printf("\nIntroduzca la aproximacion inicial:\n");
for(i=0;i scanf("\n%f", &aux);
aprox[i] = aux;}
printf("\nIntroduzca el numero de interaciones deseado:\n");
scanf("%d",&iter);
jacobi2 (matriz, v ,aprox, dimension, iter);
/* destruimos los vectores de punteros */
/*la matriz ya la libero arriba en gauss*/
free(v);
free(aprox);
system ("PAUSE");
return 0;
}

Método de Horner en Lenguaje C

#include
#include
#include

int Exponencial(float x0,int aux2){
float product;
product = 1.00;

while ((aux2)!=0){
product *= x0;
aux2--;
}
return product; /*de aqui sale el cuadrada de cada*/
}

float Horner (float x0,float *P, int N){
float *Q,aux,a;
int i,j,product;

aux = 0.00;

Q= (double*)malloc(N*sizeof(double)); /*reservamos la memoria para el vector*/


for (i=N-1;i>=0;i--){
if (i==N-1)
Q[i] = P[i]; /*bn=an*/
else
Q[i] = P[i] + (Q[i+1]*x0); /*bk=ak+ bk+1*x */

}

printf("{ ");

for ( i=0;i printf("-%f",Q[i]);
printf("} ");

j = 0;
for (i=N-1;i>=0;i--){

if (i!=0)

aux += Q[i]*Exponencial(x0,N-1-j-1);/*te calcula las b para poder multiplicarlo*/
/* a=Exponencial(x0,N-1-j-1);*/
printf("da:%f\n",aux);
j++;
}
printf ("E valor numerico del polinomio es:%f\n", aux);
return 0.00;
}
int main(){
float P[5]={-61.0,50.00,-48.00,16.00,3.00};
Horner (3,P, 5);

return 0;

}

Método de Gauss-Jacobi en Lenguaje C

#include
#include
#include


double** LeerDatos(char *nombrefichero,int *dimension)


{
double **matriz;
int i,j;
FILE *f;

if(f=fopen(nombrefichero,"r"),!f){
printf("Problema con la lectura del fichero\n");
return matriz;
}
fscanf(f,"%d\n",&(*dimension));


/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
matriz[i]=(double*)malloc((*dimension+1)*sizeof(double));


for(i=0;i<*dimension;i++)
for(j=0;j<*dimension;j++)
fscanf(f,"%lf\n",&(matriz[i][j]));


for(i=0;i<*dimension;i++)
fscanf(f,"%lf\n",&(matriz[i][*dimension]));

fclose(f);
return matriz;
}





void gauss_jacobi (double **m, float *v, float *aprox, int dimension, int iter){
int i,j,p;
float resta = 0. ,*aprox2;

aprox2 = (float*)malloc(dimension*sizeof(float));

for (i=0; i
for(p=0; p
for(j=0; j if (j!=p)
resta = resta - m[p][j] * aprox[j]; /*Se calcula el 2º termino*/

resta =resta + v[p]; /*Numerador*/
resta = resta /m[p][p];/*Denominador se divide por la diagonal(termino despejado)*/
aprox2[p]= resta;/*vector auxiliar para guardar provicionalmente hasta que termine la iteracion las aproximaciones*/
resta = 0.;
}

for (p=0;p
aprox[p]= aprox2[p]; /*En Gauss Jacobi refescamos la aproximacion al final de cada iteracion*/
}
}

for (i=0;i printf(" u[%d] = %f",i, aprox[i]);/*Mostramos las aproximaciones obtenidas con el numero de iteraciones indicado*/
}





int main(){



double **matriz;
float *v, *aprox, aux;
int dimension, iter, i, j;

matriz=LeerDatos("matriz4.dat",&dimension); /*Leemos la matriz del fichero*/

printf("\t---------------------------------------------\n");
printf("\t| METODO DE GAUSS-JACOBI |\n");
printf("\t---------------------------------------------\n");

v = (float*)malloc(dimension*sizeof(float));
aprox = (float*)malloc(dimension*sizeof(float));

printf("\n La matriz de coeficientes es: \n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",matriz[i][j]);
}
printf("}\n");
}


for (i=0;i v[i] = matriz[i][dimension];

printf("\n Introduzca la aproximacion inicial: \n");
for(i=0;i scanf("\n%f", &aux);
aprox[i] = aux;}


printf("\nIntroduzca el numero de iteraciones: \n");
scanf("%d",&iter);

printf("\n La solucion es la siguiente: \n\n");
gauss_jacobi (matriz, v ,aprox, dimension, iter);

return 0;

}

Método de Gauss-Seidel en Lenguaje C

#include
#include
#include


double** LeerDatos(char *nombrefichero,int *dimension)


{
double **matriz;
int i,j;
FILE *f;

if(f=fopen(nombrefichero,"r"),!f){
printf("Problema con la lectura del fichero\n");
return matriz;
}
fscanf(f,"%d\n",&(*dimension));


/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
matriz[i]=(double*)malloc((*dimension+1)*sizeof(double));


for(i=0;i<*dimension;i++)
for(j=0;j<*dimension;j++)
fscanf(f,"%lf\n",&(matriz[i][j]));


for(i=0;i<*dimension;i++)
fscanf(f,"%lf\n",&(matriz[i][*dimension]));

fclose(f);
return matriz;
}





void gauss_seidel (double **m, float *v, float *aprox, int dimension, int iter){
int i,j,p;
float resta = 0.;


for (i=0; i
for(p=0; p
for(j=0; j if (j!=p)
resta = resta - m[p][j] * aprox[j]; /*Se calcula el 2º termino*/

resta =resta + v[p]; /*Numerador*/
resta = resta /m[p][p];/*Denominador se divide por la diagonal(termino despejado)*/
aprox[p]= resta; /*En Gauss seidel refescamos la aproximacion en cada iteracion*/
resta = 0.; }

}

for (i=0;i printf(" u[%d] = %f",i, aprox[i]); /*Mostramos por pantalla las aproximaciones obtenidas*/
}




int main(){


double **matriz;
float *v, *aprox, aux;
int dimension, iter, i, j;

matriz=LeerDatos("matriz4.dat",&dimension); /*Leemos la matriz del fichero*/

printf("\t---------------------------------------------\n");
printf("\t| METODO DE GAUSS-SEIDEL |\n");
printf("\t---------------------------------------------\n");

v = (float*)malloc(dimension*sizeof(float));
aprox = (float*)malloc(dimension*sizeof(float));

printf("\n La matriz de coeficientes es: \n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",matriz[i][j]);
}
printf("}\n");
}


for (i=0;i v[i] = matriz[i][dimension];

printf("\n Introduzca la aproximacion inicial: \n");
for(i=0;i scanf("\n%f", &aux);
aprox[i] = aux;}


printf("\nIntroduzca el numero de iteraciones: \n");
scanf("%d",&iter);


gauss_seidel (matriz, v ,aprox, dimension, iter);
return 0;
}

Método de Gauss-Seidel en Lenguaje C (Mj_CJ)

#include
#include
#include
#define TOL 0.00001


double** LeerDatos(char *nombrefichero,int *dimension)
{
double **matriz;
int i,j;
FILE *f;

if(f=fopen(nombrefichero,"r"),!f){
printf("Problema con la lectura del fichero\n");
return matriz;
}
fscanf(f,"%d\n",&(*dimension));


/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
matriz[i]=(double*)malloc((*dimension+1)*sizeof(double));


for(i=0;i<*dimension;i++)
for(j=0;j<*dimension;j++)
fscanf(f,"%lf\n",&(matriz[i][j]));


for(i=0;i<*dimension;i++)
fscanf(f,"%lf\n",&(matriz[i][*dimension]));

fclose(f);
return matriz;
}


void liberaMatriz(double **matriz,int dimension){
int i;

for(i=0;i free(matriz[i]);
free(matriz);
}


double **multiplica(double **A,double **B,int dimension){
int i,j,k;
double **C;
C=(double**)malloc(dimension*sizeof(double*));
for(i=0;i C[i]=(double*)malloc((dimension+1)*sizeof(double));

for (i=0;i for (j=0;j for (k=0;k C[i][j]=C[i][j]+A[i][k]*B[k][j];
return C;
}

double **matrizAmpliada(double **matriz,double *vector,int dimension){
int i,j;
double **C;
C=(double**)malloc(dimension*sizeof(double*));
for(i=0;i C[i]=(double*)malloc((dimension+1)*sizeof(double));

for (i=0;i for (j=0;j C[i][j]=matriz[i][j];

for (i=0;i C[i][dimension]=vector[i];
return C;
}
double *matrizPorVector(double **matriz,double *vector,int dimension){
double *resultado,a,b,mult;
int i,j;

resultado=(double*)malloc(dimension*sizeof(double));

for(i=0;i mult=0;
for(j=0;j mult+=(matriz[i][j]*vector[j]);
resultado[i]=mult;
}

return resultado;
}


double* resolverSistemaGauss(double **A,double *b,int dimension){ /*funcion q la usa inversa*/
double M,t,*X,suma,**aux;
int i,j,p,k,c,*Fila;

Fila=(int*)malloc(dimension*sizeof(int*));
X=(double*)malloc(dimension*sizeof(double));


aux=(double**)malloc(dimension*sizeof(double*));
for(i=0;i aux[i]=(double*)malloc((dimension+1)*sizeof(double));

for(i=0;i for(j=0;j aux[i][j]=A[i][j];

aux=matrizAmpliada(aux,b,dimension);

for(j=0;j Fila[j]=j;

for(p=0;p for(k=p+1;k if(fabs(aux[Fila[k]][p])>(fabs(aux[Fila[p]][p]))){
t=Fila[p];
Fila[p]=Fila[k];
Fila[k]=t;
}

if(fabs(aux[Fila[p]][p]) < TOL){
printf("\n\n\"resolverSistemaGauss\" informa:\n");
printf("La matriz es singular\n");
return NULL;
}

for(k=p+1;k M=(aux[Fila[k]][p])/(aux[Fila[p]][p]);
for(c=p;c<=dimension;c++)
aux[Fila[k]][c]=(aux[Fila[k]][c])-(M*aux[Fila[p]][c]);
}
}

if(fabs(aux[Fila[dimension-1]][dimension-1]) < TOL){
printf("\n\n\"resolverSistemaGauss\" informa:\n");
printf("La matriz es singular\n");
return NULL;
}

X[dimension-1]=aux[Fila[dimension-1]][dimension]/(aux[Fila[dimension-1]][dimension-1]);

for(k=dimension-2;k>=0;k--){
suma=0.;
for(c=dimension-1;c>k;c--)
suma+=X[c]*aux[Fila[k]][c];
X[k]=(aux[Fila[k]][dimension]-suma)/aux[Fila[k]][k];
}
return X;
}


double **inversa(double **A,int dimension){
double **invA,*C,**aux,*b;
int i,j;

b=(double*)malloc((dimension)*sizeof(double));

invA=(double**)malloc(dimension*sizeof(double*));
for(i=0;i invA[i]=(double*)malloc((dimension)*sizeof(double));

aux=(double**)malloc(dimension*sizeof(double*));
for(i=0;i aux[i]=(double*)malloc((dimension+1)*sizeof(double));
for(i=0;i for(j=0;j aux[i][j]=A[i][j];

for(i=0;i for(j=0;j if (j==i)
b[j]=1;
else
b[j]=0;
}

if (C=resolverSistemaGauss(aux,b,dimension))
for(j=0;j invA[j][i]=C[j];
else
return NULL;
}
liberaMatriz(aux,dimension);
return invA;
}


double *suma_vector(double *A,double *B,int dimension){

int i;
double *C;

C=(double*)malloc(dimension*sizeof(double));

for (i=0;i C[i] = A[i]+B[i];

return C;
}

void Gauss_Seidel_2 (double **m, double *v, double *aprox, int dimension, int iter){
double **DL, **U, *C;
int i,j;


DL=(double**)malloc(dimension*sizeof(double*));
for(i=0;i DL[i]=(double*)malloc((dimension)*sizeof(double));

U=(double**)malloc(dimension*sizeof(double*));
for(i=0;i U[i]=(double*)malloc((dimension)*sizeof(double));


/*Hallamos la matriz D + L */

for (i=0;i for(j=0;j if (j>i)
DL[i][j] = 0.;
else
DL[i][j] = m[i][j];


/*Hallamos la inversa de D+L*/

DL = inversa(DL, dimension);

/*Obtenemos la matriz -U*/

for (i=0;i for(j=0;j if (j>i)
U[i][j] = -m[i][j];
else
U[i][j] = 0.;


/*OBTENEMOS LA MATRIZ DE GAUSS SEIDEL Mj*/

liberaMatriz(m,dimension);
m = multiplica(DL, U, dimension);
printf("\La matriz Mgs da: \n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",m[i][j]);
}
printf("}\n");
}



/*CALCULAMOS EL VECTOR Cj*/
C=matrizPorVector(DL,v,dimension);
printf("\nEl vector C da: \n");
for (j=0;j {
printf(" %4.3f ",C[j]);
}


for(i=0;i aprox = matrizPorVector(m,aprox,dimension);
aprox = suma_vector(aprox, C, dimension);}

printf("\nEl vector solucion para ese numero de iteraciones es:\n");
for(i=0;i printf("\nvector Un[%d] = %lf",i,aprox[i]);
liberaMatriz(m,dimension);
liberaMatriz(U,dimension);
liberaMatriz(DL,dimension);

}


int main(){



double **matriz;
double *v, *aprox, aux;
int dimension, iter, i, j;

matriz=LeerDatos("matriz4.dat",&dimension);

printf("\t---------------------------------------------\n");
printf("\t| METODO DE GAUSS-SEIDEL_MC |\n");
printf("\t---------------------------------------------\n");

v = (double*)malloc(dimension*sizeof(double));
aprox = (double*)malloc(dimension*sizeof(double));

printf("\n La matriz de coeficientes es: \n\n");

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",matriz[i][j]);
}
printf("}\n");
}

for (i=0;i v[i] = matriz[i][dimension];

printf("\nIntroduzca la aproximacion inicial:\n");
for(i=0;i scanf("\n%f", &aux);
aprox[i] = aux;}


printf("\nIntroduzca el numero de interaciones deseado:\n");
scanf("%d",&iter);


Gauss_Seidel_2 (matriz, v ,aprox, dimension, iter);


/* destruimos los vectores de punteros */
/*la matriz ya la libero arriba en gauss*/
free(v);
free(aprox);

return 0;

}

Método de Crout en lenguaje C

//principal.c

#include
#include
#include
#include "matrices.h"
#include "crout.h"

void captura_matriz(double **mat,int nfil,int ncol)
{
int i,j;
double valor;

for (i=0;i {
for (j=0;j {
printf("introduzca el valor[%d,%d]=",i,j);
scanf("%f",&valor);
mat[i][j]=valor;

}
}}
void muestra_matriz (double **mat,int nfil,int ncol)
{
int i,j;

for (i=0;i {
printf("{ ");
for (j=0;j {
printf("%4.3f ",mat[i][j]);
}
printf("}\n");
}
}
int main()

{
double **matriz,**matriz_u,**matriz_l,*vect,*termindp;
FILE *fich;
int i,j,tam;
/*apertura de fichero, devolviendo el tamaño de la matriz*/
tam = abrirFicheroMatriz("matriz4.dat",&fich);

/*se reserva para las matrices vectores de punteros a double que serían las filas (i)*/
matriz_u = (double**)malloc(tam*sizeof(double*));
matriz_l = (double**)malloc(tam*sizeof(double*));
vect = (double*)malloc(tam*sizeof(double));
termindp = (double*)malloc(tam*sizeof(double));
/*una vez reservadas las filas se reservan para cada puntero del vector un vector que seran
las columnas*/
for(i=0;i matriz_u[i]=(double*)malloc(tam*sizeof(double));
matriz_l[i]=(double*)malloc(tam*sizeof(double));
}
printf("\t---------------------------------------------\n");
printf("\t| METODO DE CROUT |\n");
printf("\t---------------------------------------------\n");

/*Se procedera en primer lugar a descomponer la matriz A, en matriz B y matriz B'\n");
usando para ello el metodo de Crout.\n");
Se realizara el descenso y el remonte de las misma, obtiendo asi las soluciones del sistema\n\n");*/


matriz=leerDatos(fich,tam);
/*se lee la matriz A y se muestra*/
printf("La matriz de coeficientes es: \n ");
muestra_matriz(matriz,tam,tam);

crout(matriz,matriz_u,matriz_l,tam);
/*una vez hallado B y B' o L y U son mostradas por pantalla*/

printf("\n\n El metodo de Crout concluido con exito \n");
printf("\nEsta son las dos matrices obtenidas : ");

printf("\nMatriz B' : \n ");
muestra_matriz(matriz_l,tam,tam);
printf("\nmatriz B: \n ");
muestra_matriz(matriz_u,tam,tam);
/*vector de terminos independientes*/

printf("Vector de terminos independientes :\n");
for(i=0;i {
printf("\n");
termindp[i]=matriz[i][tam];
printf("( %4.3lf )", termindp[i]);
}
printf("\n");
vect =descenso(matriz_l,termindp,tam);
vect =remonte(matriz_u,vect,tam);
printf("La resolucion del sistema de ecuaciones a concluido con exito, las soluciones son: ");
printf("\n");
for(i=0;i printf("( %4.3lf )", vect[i]);
printf("\n");
for(i=0; i {
free(matriz[i]);
free(matriz_l[i]);
free(matriz_u[i]);
}

/* destruimos losvectorres de punteros */
free(matriz);
free(matriz_l);
free(matriz_u);
free(vect);
free(termindp);
return 0;
}



//matrices.h
#ifndef MATRICES_DEF
#define MATRICES_DEF

/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A*B */
double **multiplica(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A*x, siendo x un vector */
double *matrizPorVector(double **matriz,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A+B */
double **suma(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz resultante de A-B */
double **resta(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* Le resta a la diagonal el autovalor pasado */
/* si se le pasa la matriz ampliada solo devulve la matriz */
/* sin la columna de terminos independientes */
double **matrizMenosAutovalor(double **matriz,double autovalor,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve el vector de terminos independientes de la matriz ampliada */
double *terminosIndependientes(double **matriz,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* crea una matriz ampliada pasandole una matriz y */
/* el vector de terminos independientes */
double **matrizAmpliada(double **matriz,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz traspuesta de una matriz de entrada */
double **traspuesta(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve el vector de soluciones obtenidas al realizar
el descenso en una matriz triangular superior */
double *remonte(double **B,double *b,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve el vector de soluciones obtenidas al realizar
el descenso en una matriz triangular inferior */
double *descenso(double **B,double *b,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* comprueba si una matriz es simetrica */
/* 1: SI es simetrica */
/* 0: NO es simetrica */
int simetrica(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* comprueba si dos matrices son iguales */
/* 1: SI son iguales */
/* 0: NO son iguales */
int iguales(double **A,double **B,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* comprueba si una matriz es tridiagonal */
/* 1: SI es tridiagonal */
/* 0: NO es tridiagonal */
int tridiagonal(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* imprime por pantalla la matriz A */
/* La dimension maxima para poder imprirse es 6 */
void imprimirMatrizPantalla(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* imprime a un fichero la matriz C (maximo dimension 10) */
/* requiere que el fichero haya sido abierto previamente */
void imprimirMatrizFichero(FILE *fichero,double **C,int dimension);
/*--------------------------------------------------------------------------------*/

/*--------------------------------------------------------------------------------*/
/* Analogo que imprimirMatrizPantalla pero con un vector */
void imprimirVectorPantalla(double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* Analogo que imprimirMatrizFichero pero con un vector */
void imprimirVectorFichero(FILE *fichero,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* devuelve la matriz B de la descomposicion A=B*Bt */
/* comprueba que la matriz A sea simetrica, devolviendo
un NULL en caso de no ser simetrica */
double **descomposicionCholesky (double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* lee de un fichero de texto una matriz */
/* la matriz tiene que venir de la siguiente manera */
/* en la primera linea la dimension, despues en cada linea
se van introduciendo los elementos de las distintas filas
ordenadamente. Al final y tambien en cada linea los terminos
independientes; Ejemplo:
2
1
0
0
1
1
1
Esta seria la matriz identidad de orden dos y con el vector (1,1) como
terminos independientes. */
double** leerDatos(FILE *f,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* abre un fichero que contiene una matriz con la representacion
anterior devolviendo:
la dimension si pudo abrir el fichero satisfactoriamente
0 en caso contrario */
int abrirFicheroMatriz(char *nombreFichero,FILE **fentrada);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* crea un fichero con una matriz con la representacion comentada en leerDatos
para luego poder se leida por leerDatos. Como parametros requiere la matriz,
el vector de terminos independientes y la dimension de la matriz.
Devuelve:
1: pudo crear el fichero correctamente
0: no pudo crear el fichero */
int crearFicheroMatriz(double **matriz,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* elimina el espacio reservado para una matriz */
void eliminaMatriz(double **matriz,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* Esta funcion se apoya en la funcion descomposicionCholesky para obtener la
/* matriz B y luego resolver el sistema B*z=b mediante un remonte. La solucion
/* de este sistema la utilizamos para resolver el sistema Bt*u=z mediante un
/* descenso. Requiere que la matriz sea definida positiva y simetrica */
double *resolverSistemaCholesky(double **A,double *vector,int dimension);
/*--------------------------------------------------------------------------------*/


double* resolverSistemaGauss(double **A,double *b,int dimension);

double* resolverSistemaGaussSeidel(double **A,double *b,int dimension,int *Nmax);

double *autovaloresJacobi(double **matriz,int dimension,int *Nmax);

double determinante(double **matriz,int dimension);


/*--------------------------------------------------------------------------------*/
/* Calcula la inversa de una matriz. Para ello resuelve tantos sistemas como sea
/* la dimension de entrada. La peculiaridad de estos sistemas es que tiene como
/* terminos independientes las columnas de la matriz identidad. Los vectores obtenidos
/* como soluciones de estos sistemas seran las columnas de la matriz inversa. Para
/* resolver los distintos sistemas utilizamos la funcion resolverSistemaGauss. Si
/* Gauss no es capaz de resolver el sistema devuelve NULL */
double **inversa(double **A,int dimension);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* calcula el autovalor maximo de la matriz de entrada que en valor absoluto es */
/* el radio espectral de la matriz de entrada */
double potencia(double **matriz, double *x,double *autovector,int dimension,int *Nmax);
/*--------------------------------------------------------------------------------*/


/*--------------------------------------------------------------------------------*/
/* resuelve un sistema de ecuaciones por el metodo de la relajacion, que es muy
/* parecido al metodo de Gauss-Seidel pero introduciendo un parametro de
/* relajacion w. Para el calculo de este w nos apoyamos en el metodo de la
/* potencia para el calculo del autovalor maximo. Para ello es imprescindible
/* que la matriz de entrada sea tridiagonal. En caso de no ser tridiagonal
/* resolvemos el sistema con una w predefinida (ver implementacion) */
double* resolverSistemaRelajacion(double **A,double *b,int dimension,int *Nmax);
/*--------------------------------------------------------------------------------*/

double potenciaInversa(double **matriz, double *uInicial,double *autovector,int dimension,int *Nmax);

#endif

//matrices.c

#include
#include
#include
#include
#include "matrices.h"
#define TOL 0.00000001

int iguales(double **A,double **B,int dimension){
int i,j;
for(i=0;i for(j=0;j if(fabs((A[i][j])-(B[i][j])) > TOL)
return 0;
return 1;
printf("\nLas matrices estan bien");
}

int abrirFicheroMatriz(char *nombreFichero,FILE **fentrada){
int dimension;
if(*fentrada=fopen(nombreFichero,"rb"),!*fentrada){
printf("\n\n\"abrirFicheroMatriz\" informa:\n");
printf("ERROR: No se ha podido abrir el fichero \"%s\"",nombreFichero);
return 0;
}
fscanf(*fentrada,"%d\n",&dimension);
return dimension;
}

double **multiplica(double **A,double **B,int dimension){
int i,j,k;
double **C;
C=(double**)malloc(dimension*sizeof(double*));
for(i=0;i C[i]=(double*)malloc((dimension+1)*sizeof(double));
for (i=0;i for (j=0;j for (k=0;k C[i][j]=C[i][j]+A[i][k]*B[k][j];
return C;
}

double** leerDatos(FILE *f,int dimension){
double **matriz;
int i,j;

/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(dimension*sizeof(double*));
for(i=0;i matriz[i]=(double*)malloc((dimension+1)*sizeof(double));

for(i=0;i for(j=0;j fscanf(f,"%lf\n",&(matriz[i][j]));

for(i=0;i fscanf(f,"%lf\n",&(matriz[i][dimension]));

fclose(f);
return matriz;
}

int tridiagonal(double **A,int dimension){
int i,j;
for(i=2;i for(j=0;j if((fabs(A[j][i])+fabs(A[i][j])) > TOL)
return 0;
return 1;
}

double *descenso(double **B,double *b,int dimension){
int i,j;
double suma,*vec;
vec=(double*)malloc(dimension*sizeof(double));
vec[0]=b[0]/(B[0][0]);
for(i=1;i suma=0.;
for(j=0;j suma+=vec[j]*B[i][j];
vec[i]=(b[i]-suma)/B[i][i];
}
return vec;
}

double *remonte(double **B,double *b,int dimension){
int i,j;
double suma,*vec;
vec=(double*)malloc(dimension*sizeof(double));
vec[dimension-1]=b[dimension-1]/(B[dimension-1][dimension-1]);
for(i=dimension-2;i>=0;i--){
suma=0.;
for(j=dimension-1;j>i;j--)
suma+=vec[j]*B[i][j];
vec[i]=(b[i]-suma)/B[i][i];
}
return vec;
}

//crout.h

#ifndef CROUT
#define CROUT

/**************METODO DE FACTORIZACION DE CROUT*****************/

void crout(double **mat,double **matriz_u,double **matriz_l,int tam);
/*procedimiento que dada una matriz y un tamaño devuelve dos matrices triangulares que su multiplicacion dan mat, dado que crout es un metodo de factorizacion*/
#endif


//crout.c

void crout(double **mat,double **mat_u,double **mat_l,int tam)
{

int i,j,k;


for (i=0; i {
for(j=0;j {
if (i==j)
mat_u[i][j]=1.;
else
mat_u[i][j]=0.;

mat_l[i][j]=0.;
}
}

mat_l[0][0]=mat[0][0];
mat_u[0][1]=(mat[0][1]/mat_l[0][0]);
mat_l[1][0]=mat[1][0];

/*Se rellenan las matrices*/

for(i=1;i {
mat_l[i][i]=mat[i][i]-(mat_l[i][i-1]*mat_u[i-1][i]);
mat_u[i][i+1]=(mat[i][i+1]/mat_l[i][i]);
mat_l[i+1][i]=mat[i+1][i];
}
mat_l[tam-1][tam-2]=mat[tam-1][tam-2];
mat_l[tam-1][tam-1]=mat[tam-1][tam-1]-mat_l[tam-1][tam-2]*mat_u[tam-2][tam-1];
}

Método de la Bisección en Lenguaje C

#include
#include

#define Cero 0.00001

float f(float x){
return (x*x*x)-27;
}

float biseccion (float a,float b, float tol){
float c = (a+b)/2.0;

while ( (f(c) != Cero) && ( fabs(b-a) > tol ) ){
if (f(a) * f(c) < Cero)
b=c;
else
a=c;
c = (a+b)/2.0;}
return c;
}

int main(){

printf("el valor de la raiz es: %f \n", biseccion (0.0, 3.0, 0.001) );
return 0;

}

Método de Cholesky en Lenguaje C

//Principal.c
//Autores.- Sergio Romero Sánchez & Alejandro E. Castro Robredo
//Práctica Método de Cholesky.-
//Especialidad.- I.T.I.G.
//Fecha.- 18/01/2005

#include
#include
#include "cholesky.h"
#include "Leer.h"
int main()
{
double **A,*x;
int n,i,j;

A=LeerDatos("Dto2.dat",&n);
if ( (x=Cholesky(A,n))!=NULL)
{
printf ("\nEstos son los resultados: \n");
for(i=0;i printf("\n X%d = %.14lf",i,x[i]);
free(x);
printf ("\n");
}
else
{
printf ("No se ha podido resolver el sistema \n");
}
for (i=0;i {
free(A[i]);
}
free(A);
return 1;
}




//Leer.h
double** LeerDatos(char *nombrefichero,int *dimension);

//Leer.c

#include
#include "Leer.h"
#include

double** LeerDatos(char *nombrefichero,int *dimension) {
double **matriz;
int i,j;
FILE *f;

if(f=fopen(nombrefichero,"r"),!f){
printf("Problema con la lectura del fichero\n");
return matriz;
}
fscanf(f,"%d\n",&(*dimension));

/* Asignacion de memoria dinamica*/
matriz=(double**)malloc(*dimension*sizeof(double*));
for(i=0;i<*dimension;i++)
matriz[i]=(double*)malloc((*dimension+1)*sizeof(double));

for(i=0;i<*dimension;i++)
for(j=0;j<*dimension;j++)
fscanf(f,"%lf\n",&(matriz[i][j]));

for(i=0;i<*dimension;i++)
fscanf(f,"%lf\n",&(matriz[i][*dimension]));

fclose(f);
return matriz;
}



//Cholesky.h

int Diagonalizar(double **A,double **B,int dimension);
void MatrizTraspuesta(double **A,double **At,int dimension);
void DescensoCholesky(double **B,double *y,double *b,int dimension);
void AscensoCholesky(double **Bt,double *x,double *y,int dimension);
int Simetrica (double **A,int n);
double** LeerDatos(char *nombrefichero,int *dimension);
double* Cholesky (double **A, int dimension);

//Cholesky.c

#include "cholesky.h"
#include
#include
#include

int Diagonalizar(double **A,double **B,int dimension)
/* Ambas matrices deben ya disponer de espacio de memoria dinamica
// A = matriz original de coeficientes
// B = triangular inferior (los de arriba todos a cero) despues de
aplicar Cholesky
// dimension = dimension de la matriz cuadrada*/
{
int i,j,k;
double suma,producto,z;

for(j=0;j {
for (i=j;i {
if(i==j)
{
suma = 0.;
for(k=0;k {
producto = B[j][k]*B[j][k];
suma = suma + producto;
}
if ((z=A[j][j]-suma)<=0)
{
return 0;
}
else
{
B[j][j] = sqrt(z);
}
}
else
{
suma = 0.;
for(k=0;k {
producto = B[i][k]*B[j][k];
suma = suma + producto;
}
B[i][j] = (double)(A[i][j]-suma)/B[j][j];
}
printf("\nB[%d][%d] = %.14lf",i,j,B[i][j]);
}
}
return 1;
}

void MatrizTraspuesta(double **A,double **At,int dimension)
/* Ambas matrices deben ya disponer de espacio de memoria dinamica
// A = matriz original
// At = matriz traspuesta de la matriz A
// dimension = dimension de la matriz cuadrada*/
{
int i,j;

for (j=0;j {
for(i=0;i At[i][j] = A[j][i];
}
return;
}


void DescensoCholesky(double **B,double *y,double *b,int dimension)
/*Ambas matrices deben ya disponer de espacio de memoria dinamica
// B = matriz B
// y = soluciones del descenso
// b = vector de terminos independientes
// dimension = dimension de la matriz cuadrada*/
{
int i,j;
double suma;

for(i=0;i {
suma = 0.;
for(j=0;j {
suma = suma + (B[i][j]*y[j]);
}
y[i] = (double)(b[i]-suma)/B[i][i];
}

for(i=0;i<4;i++)
printf("\ny%d = %.14lf",i,y[i]);

return;
}

void AscensoCholesky(double **Bt,double *x,double *y,int dimension)
/*Ambas matrices deben ya disponer de espacio de memoria dinamica
// Bt = matriz traspuesta de B
// x = soluciones del sistemas
// y = soluciones del descenso
// dimension = dimension de la matriz cuadrada */
{
int i,j;
double suma;

for(i=(dimension-1);i>=0;i--)
{
suma = 0.;
for(j=i;j {
suma = suma + (Bt[i][j]*x[j]);
}
x[i] = (double)(y[i]-suma)/Bt[i][i];
}
return;
}
int Simetrica (double **A,int n)

/* A=Matriz de coeficientes si es simétrica devuelve 1, en caso contrario
devuelve 0 */
{
int i,j;

for (i=0;i {
for (j=i+1;j {
if (A[i][j]!=A[j][i])
return 0;
}
}
return 1;
}

double* Cholesky (double **A, int dimension)
{
/* Matriz A contiene coefcientes y tréminos independientes
dimension de la matriz
Si A cumple condiciones de Cholesky entonces
Devuelve X con la solucion del sistema
en caso contrario
Devuelve NULL */

double *b,**B,**Bt,*x,*y;
int h,i;
if (Simetrica(A,dimension)==1)
{
b=(double*)malloc(dimension*sizeof(double));
B = (double**)malloc(dimension*sizeof(double*));
for(h=0;h B[h] = (double*)malloc(dimension*sizeof(double));

Bt = (double**)malloc(dimension*sizeof(double*));
for(h=0;h Bt[h] = (double*)malloc(dimension*sizeof(double));
x = (double*)malloc(dimension*sizeof(double));
y = (double*)malloc(dimension*sizeof(double));
for(i=0;i b[i]=A[i][dimension];

if (Diagonalizar(A,B,dimension)==1)
{
MatrizTraspuesta(B,Bt,4);
DescensoCholesky(B,y,b,dimension);
AscensoCholesky(Bt,x,y,dimension);

for (i=0;i {
free(B[i]);free(Bt[i]);
}
free(y);
free(B);
free(Bt);
free(b);
return x;
}
else
{
printf ("La matriz de coeficientes no es definida positiva,y por
tanto no se puede resolver a través de Cholesky \n");
for (i=0;i {
free(B[i]);free(Bt[i]);
}
free(y);free(B);
free(Bt);free(b);
return NULL;
}}
else
{
printf ("La matriz de coeficientes no es simétrica,y por tanto no
se puede resolver a través de Cholesky \n");
return NULL;
}}