[C] Integrazione di Gauss Legendre.

Suppish · Messaggio da **Suppish** » martedì 22 aprile 2014, 14:37

Salve,
Devo creare un programma che mi calcoli l'integrale di una certa funzione con lo sviluppo di Gauss-Legendre.
La teoria che c'è dietro (in pochissime parole) afferma che l'integrale della funzione data si può scrivere come Immagine

Il programma è notevolmente semplificato giacchè ho a disposizione le coordinate dei punti ξ e i pesi w su file. Resta solo da trovare un modo efficiente di chiamare i giusti valori di pesi e punti tra quelli disponibili, in base al numero n di punti su cui si vuole estendere quella sommatoria.
Il file dei pesi e dei punti è il seguente (per semplicità è richiesto che contempli i casi fino a n = 16)
http://www.holoborodko.com/pavel/numeri ... tegration/

In pratica, se l'utente ad esempio sceglie n = 4, vorrei che il programma leggesse solo i valori per n = 4 ignorando gli altri e che si adattasse per eseguire la sommatoria a queste condizioni (quindi ciclare per i che va da 0 a 4 i vari valori di punti e pesi letti dal file).
Per ora l'unica soluzione da me trovata è quella di scrivere tutti i casi possibili e in base al valori di punti scelto dall'utente andare a eseguire esclusivamente il caso specificato.
Come potete vedere è una soluzione grezza e lenta. Richiederebbe pesanti modifiche laddove volessi contemplare casi per n>16.
Ci sono altre maniere?

Grazie per l'attenzione,
Fabio

ienaplinsky

Ciao scusa ma su quel file ad alcuni valori sono associate più ascisse bisogna prenderle tutte ?
Che intendi per casi possibili?
Posta quello che hai fatto

Suppish · Messaggio da **Suppish** » martedì 22 aprile 2014, 17:17

Si, il programma deve eseguire il conto in formula per tutte le ascisse e i rispettivi pesi per quel dato n.

Casi possibili intendo n possibili. La scelta di n è dell'utente.

Fra qualche ora, appena riprendo possesso del pc fisso a casa, posto anche una bozza del programma nella soluzione che avevo pensato.

ienaplinsky

servirebbe capire anche come è fatto il file cioè come hai impostato il file o se hai dei limiti nel impostazione del file in modo da capire come andare a leggere al meglio i dati
se puoi ad esempio impostare i dati nel file in questo modo
la prima colonna è n la seconda l' ascissa e la terza il peso (non so se ci vuole sia + che meno ho preso solo i positivi ma non cambia molto)

Codice: Seleziona tutto

2 0.5773502691896257645091488 1.0000000000000000000000000
3 0 0.8888888888888888888888889
3 0.7745966692414833770358531 0.5555555555555555555555556
4 0.3399810435848562648026658 0.6521451548625461426269361
4 0.8611363115940525752239465 0.3478548451374538573730639

puoi leggere velocemente e mentre leggi fare anche i calcoli

Codice: Seleziona tutto

int main(void) {
    FILE *file = fopen("file", "r");
    int i, n;
    double ascissa, peso;
    
    if(!file) {
        puts("Errore nell' apertura del file");
        return -1;
    }
    
    do {
        puts("inserire n");
        scanf("%d", &n);
    } while(n < 2 || n > 16);
    
    while(fscanf(file, "%d", &i) != EOF && i <= n) {
        fscanf(file, "%lf", &ascissa);
        fscanf(file, "%lf", &peso);
        // fai i conti
    }
    
    return 0;
}

ps: vedi che nei double non credo ci entrino numeri cosi precisi

Vincenzo1968

Ho provato la libreria del sito di Holoborodko:

Sembra veloce:

Codice: Seleziona tutto

[vincenzo]$ gcc -Wall -W -pedantic -O2 gauss_legendre.c example.c -o example -lm
                          
[vincenzo]$ time ./example
Numerical Approximation of int(sin(x), x=0..Pi) by Gauss-Legendre Quadrature:
n =    2:	approx value = 1.935819574651137298105	error = 0.0641804253488627
n =    3:	approx value = 2.001388913607743180734	error = 0.00138891360774318
n =    4:	approx value = 1.999984228457722057470	error = 1.57715422779425e-05
n =    5:	approx value = 2.000000110284471777078	error = 1.10284471777078e-07
n =    6:	approx value = 1.999999999477270362647	error = 5.22729637353336e-10
n =    7:	approx value = 2.000000000001790567694	error = 1.79056769411545e-12
n =    8:	approx value = 1.999999999999995115019	error = 4.88498130835069e-15
n =    9:	approx value = 2.000000000000000000000	error = 0
n =   10:	approx value = 2.000000000000000000000	error = 0
n =   11:	approx value = 2.000000000000000000000	error = 0
n =   12:	approx value = 2.000000000000000000000	error = 0
n =   13:	approx value = 2.000000000000000000000	error = 0
n =   14:	approx value = 2.000000000000000000000	error = 0
n =   15:	approx value = 2.000000000000000000000	error = 0
n =   16:	approx value = 2.000000000000000000000	error = 0
n =   17:	approx value = 2.000000000000000444089	error = 4.44089209850063e-16
n =   18:	approx value = 2.000000000000000444089	error = 4.44089209850063e-16
n =   19:	approx value = 2.000000000000000000000	error = 0
n =   20:	approx value = 2.000000000000000444089	error = 4.44089209850063e-16
n =   21:	approx value = 1.999999999999999333866	error = 6.66133814775094e-16
n =   22:	approx value = 1.999999999999999777955	error = 2.22044604925031e-16
n =   23:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   24:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   25:	approx value = 1.999999999999999777955	error = 2.22044604925031e-16
n =   26:	approx value = 1.999999999999999777955	error = 2.22044604925031e-16
n =   27:	approx value = 1.999999999999999333866	error = 6.66133814775094e-16
n =   28:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   29:	approx value = 1.999999999999999777955	error = 2.22044604925031e-16
n =   30:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   31:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   32:	approx value = 1.999999999999999777955	error = 2.22044604925031e-16
n =   33:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   34:	approx value = 1.999999999999999333866	error = 6.66133814775094e-16
n =   35:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   36:	approx value = 2.000000000000000000000	error = 0
n =   37:	approx value = 1.999999999999999333866	error = 6.66133814775094e-16
n =   38:	approx value = 1.999999999972282616056	error = 2.77173839435818e-11
n =   39:	approx value = 1.999999999999999777955	error = 2.22044604925031e-16
n =   40:	approx value = 1.999999999999998667732	error = 1.33226762955019e-15
n =   41:	approx value = 2.000000000019693136011	error = 1.9693136010801e-11
n =   42:	approx value = 1.999999999999997335465	error = 2.66453525910038e-15
n =   43:	approx value = 1.999999999986464604973	error = 1.35353950270201e-11
n =   44:	approx value = 1.999999999996994182183	error = 3.00581781687015e-12
n =   45:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   46:	approx value = 2.000000000006650680007	error = 6.65068000671454e-12
n =   47:	approx value = 1.999999999999998889777	error = 1.11022302462516e-15
n =   48:	approx value = 1.999999999991228127882	error = 8.77187211756336e-12
n =   49:	approx value = 2.000000000019829471398	error = 1.9829471398225e-11
n =   50:	approx value = 1.999999999984244603013	error = 1.57553969870605e-11
n =   51:	approx value = 2.000000000000223376873	error = 2.23376872554581e-13
n =   52:	approx value = 2.000000000006295852728	error = 6.29585272804434e-12
n =   53:	approx value = 1.999999999992864152532	error = 7.13584746847573e-12
n =   54:	approx value = 1.999999999986300291965	error = 1.36997080346646e-11
n =   55:	approx value = 1.999999999989424459557	error = 1.05755404433694e-11
n =   56:	approx value = 1.999999999980550668965	error = 1.94493310345933e-11
n =   57:	approx value = 1.999999999994776622714	error = 5.22337728625644e-12
n =   58:	approx value = 1.999999999969677588751	error = 3.03224112485623e-11
n =   59:	approx value = 1.999999999972649877833	error = 2.73501221670358e-11
n =   60:	approx value = 1.999999999986062926283	error = 1.39370737173294e-11
n =   61:	approx value = 1.999999999951697526868	error = 4.83024731323667e-11
n =   62:	approx value = 1.999999999960391461329	error = 3.9608538671132e-11
n =   63:	approx value = 1.999999999971250552733	error = 2.87494472672734e-11
n =   64:	approx value = 2.000000000000000000000	error = 0
n =   65:	approx value = 1.999999999970391906245	error = 2.96080937545184e-11
n =   66:	approx value = 1.999999999979707343556	error = 2.02926564440986e-11
n =   67:	approx value = 1.999999999990506704961	error = 9.49329503896479e-12
n =   68:	approx value = 1.999999999977978060173	error = 2.20219398272548e-11
n =   69:	approx value = 1.999999999985965892790	error = 1.40341072096817e-11
n =   70:	approx value = 1.999999999995084598581	error = 4.91540141922542e-12
n =   71:	approx value = 1.999999999983752108079	error = 1.62478919207842e-11
n =   72:	approx value = 1.999999999990609955702	error = 9.39004429767465e-12
n =   73:	approx value = 1.999999999998333555240	error = 1.66644475996236e-12
n =   74:	approx value = 2.000000000006834532940	error = 6.83453293959246e-12
n =   75:	approx value = 1.999999999994054977748	error = 5.94502225226279e-12
n =   76:	approx value = 2.000000000000623945340	error = 6.23945339839338e-13
n =   77:	approx value = 2.000000000007795097900	error = 7.79509790049815e-12
n =   78:	approx value = 2.000000000015492940264	error = 1.54929402640391e-11
n =   79:	approx value = 2.000000000002217781514	error = 2.21778151399121e-12
n =   80:	approx value = 2.000000000008296474618	error = 8.29647461841887e-12
n =   81:	approx value = 2.000000000014788170688	error = 1.47881706880071e-11
n =   82:	approx value = 2.000000000021638690839	error = 2.16386908391542e-11
n =   83:	approx value = 2.000000000008482992087	error = 8.4829920865559e-12
n =   84:	approx value = 2.000000000013987477843	error = 1.39874778426474e-11
n =   85:	approx value = 2.000000000019776180693	error = 1.9776180693043e-11
n =   86:	approx value = 2.000000000025802471271	error = 2.58024712707083e-11
n =   87:	approx value = 2.000000000013144152433	error = 1.31441524331422e-11
n =   88:	approx value = 2.000000000018061108165	error = 1.8061108164602e-11
n =   89:	approx value = 2.000000000023169466346	error = 2.31694663455073e-11
n =   90:	approx value = 2.000000000028435476196	error = 2.84354761959094e-11
n =   91:	approx value = 2.000000000033829383739	error = 3.38293837387482e-11
n =   92:	approx value = 2.000000000020845991600	error = 2.08459915995718e-11
n =   93:	approx value = 2.000000000025327295816	error = 2.53272958161688e-11
n =   94:	approx value = 2.000000000029910740551	error = 2.99107405510313e-11
n =   95:	approx value = 2.000000000034575453611	error = 3.45754536112963e-11
n =   96:	approx value = 2.000000000000000444089	error = 4.44089209850063e-16
n =   97:	approx value = 2.000000000026541879805	error = 2.65418798051087e-11
n =   98:	approx value = 2.000000000030523583661	error = 3.05235836606244e-11
n =   99:	approx value = 2.000000000034553693240	error = 3.45536932400137e-11
n =  100:	approx value = 1.999999999999998667732	error = 1.33226762955019e-15
n =  101:	approx value = 2.000000000042685410762	error = 4.26854107615782e-11
n =  102:	approx value = 2.000000000046756376548	error = 4.67563765482737e-11
n =  103:	approx value = 2.000000000033986591319	error = 3.39865913190351e-11
n =  104:	approx value = 2.000000000037481573401	error = 3.74815734005551e-11
n =  105:	approx value = 2.000000000040976111393	error = 4.09761113928653e-11
n =  106:	approx value = 2.000000000044461323512	error = 4.44613235117686e-11
n =  107:	approx value = 2.000000000047922110724	error = 4.79221107241301e-11
n =  108:	approx value = 2.000000000051355808495	error = 5.13558084946908e-11
n =  109:	approx value = 2.000000000039064307344	error = 3.90643073444608e-11
n =  110:	approx value = 2.000000000042054359994	error = 4.20543599943812e-11
n =  111:	approx value = 2.000000000045019099559	error = 4.50190995593402e-11
n =  112:	approx value = 2.000000000047952308790	error = 4.79523087903999e-11
n =  113:	approx value = 2.000000000050846882260	error = 5.08468822602026e-11
n =  114:	approx value = 2.000000000053698379077	error = 5.36983790766499e-11
n =  115:	approx value = 2.000000000056498805634	error = 5.64988056339644e-11
n =  116:	approx value = 2.000000000059248161932	error = 5.92481619321461e-11
n =  117:	approx value = 2.000000000061937566187	error = 6.19375661869981e-11
n =  118:	approx value = 2.000000000049603432473	error = 4.96034324726224e-11
n =  119:	approx value = 2.000000000051995741046	error = 5.19957410460847e-11
n =  120:	approx value = 2.000000000054347637501	error = 5.43476375014507e-11
n =  121:	approx value = 2.000000000056644910984	error = 5.6644910984005e-11
n =  122:	approx value = 2.000000000058886229226	error = 5.88862292261183e-11
n =  123:	approx value = 2.000000000061074700852	error = 6.10747008522594e-11
n =  124:	approx value = 2.000000000063203664524	error = 6.32036645242806e-11
n =  125:	approx value = 2.000000000065271787975	error = 6.52717879745524e-11
n =  126:	approx value = 2.000000000067279515292	error = 6.72795152922845e-11
n =  127:	approx value = 2.000000000069224626031	error = 6.92246260314278e-11
n =  128:	approx value = 1.999999999999999777955	error = 2.22044604925031e-16
--------------------------------------------------------------------------------
n =  511:	approx value = 2.000000000017685852782	error = 1.76858527822787e-11
n =  512:	approx value = 1.999999999999999333866	error = 6.66133814775094e-16
n =  513:	approx value = 2.000000000017450041412	error = 1.74500414118484e-11

real	0m0.143s
user	0m0.140s
sys	0m0.000s

Ho modificato leggermente il file example.c:

Codice: Seleziona tutto

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

#ifndef PI
	#define PI 3.1415926535897932384626433832795028841971693993751
#endif
#ifndef FABS
	#define FABS(a) ((a)>=0?(a):-(a))
#endif

double f(double x, void* data)
{
	return sin(x);
}

int main( /*int argc, char* argv[] */)
{

	/* numerical approximation of integral */
	double approx;		

	/* true value of int(sin(x), x=0..Pi) = 2.0*/
	double exact = 2.0; 

	/* approximation error */
	double error;       

	int i;

	printf("Numerical Approximation of int(sin(x), x=0..Pi) by Gauss-Legendre Quadrature:\n");
	for (i = 2; i <= 128; i++)
	{
		approx = gauss_legendre(i, f, NULL, 0, PI);
		error = approx - exact;
		/* printf("n = %4d: error = %.15g\n",i,FABS(error)); */
		printf("n = %4d:\tapprox value = %15.21f\terror = %.15g\n", i, approx, FABS(error));
	}
	
	printf("--------------------------------------------------------------------------------\n");

	i = 511;
	approx = gauss_legendre(i, f, NULL, 0, PI);
	error = approx - exact;
	printf("n = %4d:\tapprox value = %15.21f\terror = %.15g\n", i, approx, FABS(error));

	i = 512;
	approx = gauss_legendre(i, f, NULL, 0, PI);
	error = approx - exact;
	printf("n = %4d:\tapprox value = %15.21f\terror = %.15g\n", i, approx, FABS(error));
	
	i = 513;
	approx = gauss_legendre(i, f, NULL, 0, PI);
	error = approx - exact;
	printf("n = %4d:\tapprox value = %15.21f\terror = %.15g\n", i, approx, FABS(error));
	
		
	return 0;
}

Forum Ubuntu-it