🚀 go-pugleaf

-----BEGIN PGP SIGNED MESSAGE-----


    >>>> El Wed, 5 Apr 2000 23:58:35 +0200, un alien�gena
    >>>> intent� hacerse pasar por gradha@iname.com,
    >>>> y escribi� en es.comp.os.linux.programacion:

>Saludos.
>
>Es posible que dentro de poco vaya a tener que hacer un proyecto para la
>universidad, de estad�stica, c�lculo de integrales y dem�s cosas "raras" :-)
>
>El caso es que a pesar de haber nacido con un puntero bajo el brazo no
>tengo por d�nde cojer el hacer una integral en C, o cualquier otra cosa
>de ese estilo, porque nunca lo he intentado antes y porque soy tremendamente
>malo en asignaturas claramente no-inform�ticas.
>
>Ya que este no es el grupo, �d�nde puedo conseguir informaci�n sobre c�mo
>atacar el tema? Preferir�a que fuese sin c�digo, a nivel te�rico, tipo
>qu� acciones tomar ante una integral de no-se-qu� con tal par�metro o
>documentos que traten sobre la p�rdida de precisi�n en los decimales,
>si hay alguna forma inform�tica o matem�tica para evitarla, etc.

Puedes usar una librer�a de aritm�tica multiprecisi�n como GMP y olvidarte de
los detalles num�ricos.

Por si te interesa aqu� va un algoritmo de integraci�n adaptativa, mediante
f�rmula de cuadratura o integraci�n v�a interpolaci�n polinomial, usando la
regla de Simpson.

Jeje, siempre es �til tener a mano la chuleta... bueno, aqu� va el
pseudo-c�digo (si quieres una explicaci�n y justificaci�n matem�tica te la doy
otro d�a):

Antes de nada, definimos la funci�n en f (en la pr�ctica ser� pasada como
par�metro a la funci�n de integraci�n, como puntero a funci�n), a es el
comienzo del intervalo, b el extremo derecho, e una acotaci�n del error, y n
el n�mero m�ximo de subdivisiones del intervalo.

input a,b,e,n
begin
	delta=b-a; suma=0; h=delta/2; c=(a+b)/2; k=1;
	fa=f(a); fb=f(b); fc=f(c);
	S=(fa+4*fc+fb)*h/3;
	v[k] = [a,h,fa,fc,fb,S];
	while 1 <= k <= n
		h=v[k][2]/2;
		fy=f(v[k][1]+h);
		S1=(v[k][3]+4*fy+v[k][4])*h/3;
		fz=f(v[k][1]+3*h);
		S2=(v[k][4]+4*fz+v[k][5])*h/3;
		if abs(S1+S2-v[k][6]) < 30*e*h/delta then
			suma=suma+S1+S2+1/15*(S1+S2-v[k][6]);
			k=k-1;
			if k<=0 then output suma; stop;
		else
			fv=v[k][5];
			v[k]= [v[k][1],h,v[k][3],fy,v[k][4],S1];
			k=k+1;
			v[k]=[v[k-1][1]+2*h,h,v[k-1][5],fz,fv,S2];
		endif
	endwhile
	output fallo
end


donde v[k] se usa a modo de pila, cada elemento tiene 6 componentes:
	v[k] = [ v[k][1],v[k][2],.... ]
Esto es lo que significa cada componente:
1 punto a la izquierda del intervalo
2 mitad de la longitud del intervalo
3 valor de f en el extremo izquierdo
4    "             punto medio    
5    "             extremo derecho
6 valor que la regla de Simpson da a este intervalo

Lo que el pseudoc�digo anterior hace es calcular el valor aproximado de la
integral en el segmento mediante la regla de Simpson: 
         int ( f(x), x=a..b ) ~= (b-a)/6*(f(a)+4*f((a+b)/2)+f(b))

despu�s dividimos el intervalo en dos subintervalos, parti�ndolo por la mitad,
y sacamos el valor estimado por Simpson de las dos integrales
correspondientes. El error se estima comparando la primera aproximaci�n con
la suma de las dos aproximaciones de los subintervalos. Si es menor que e,
hemos acabado.

- -- 
              _           _                             
             | |__   __ _| |_ ___ _ __ ___   __ _ _ __  
             | '_ \ / _` | __/ __| '_ ` _ \ / _` | '_ \ 
             | |_) | (_| | |_\__ \ | | | | | (_| | | | |
             |_.__/ \__,_|\__|___/_| |_| |_|\__,_|_| |_|

             Real email: batsman dot geo at yahoo dot com

-----BEGIN PGP SIGNATURE-----
Version: 2.6.3ia
Charset: noconv
Comment: Public key via email: send mail with Subject: [PGPKEY]

iQEVAwUBOOz3QL0Y27yr0iefAQHNOgf/WyX7VwGyHCHIn7s3JgS53GOamjr4uFF1
Uz1mWmLqyYyQhEH/nNLxpTxnD9oFhNCKqs2bIhWd0bE1XTJ0mG46G/iB66NuGvnL
bVFzFqf2hFWytAmpywjk/5t63O9HgIs+oev2VJ9zExz2ZwAik0QqeD1aQNB9QkEA
xqPZrOianpWqTteVIjzO+LMY9E5mROngHHRuKqKTEcwl+oqDVe7aLbeDCEgydAjD
JL1mQtErgxJDe/rOMXn34MTXLXcYIZqRaRfA0ArVSILP+fJB8JQExKZtNBPoo0Ze
J7f52scvDezba1A7OzswjyR0A3eH0AxwQQXYW/2GOPQNqvgwrLZU/w==
=bV+R
-----END PGP SIGNATURE-----