Calcolo sommatoria

[irelimax]

Ciao a tutti,

sto implementando un programma che mi calcola la sommatoria della formula di Bessel:
\(\displaystyle \sum_{k=0}^{ \infty}{\frac{(-1)^k(x/2)^{2k}}{k!(n+k)!}}) \)

Il problema è che facendomi stampare tutte le somme successive, si vede che la sommatoria a partire dal 5°/6° step mi restituisce sempre lo stesso numero. Immagino sia un problema di casting ma non sò come risolverlo.

Qualcuno potrebbe aiutarmi?


Il programma che faccio girare è il seguente:

#include <conio.h>
#include <math.h>
#include <iostream>
#include<fstream>
using namespace std;
 
double fact (int n) //factorial function
{
	double F=1.0;
	for (int i=1; i<=(double)n; i++)
	{
		F=F*i;
	}
	return F;
}



double Sum (int n, double x) //evaluation Bessel sum
{
	double sgn=1.0, num=1.0, den=fact(n);
	double S=1.0/(fact(n));
	int N=100;
	for(int i=1; i<N; i++)
	{
	sgn=-sgn;
	num*=(x/2.0)*(x/2.0);
    den*=(double)i*(n+i);
	S+=sgn*(num/den);
	cout<<"S"<<i<<"="<<S<<endl;
	}

	S+=pow(x/2,(double) n);
	return S;
}



int main()
{
	int n;
	double x,S;
	cout<<"Please, enter the order n and the value of x to evaluate numerically the Bessel function J_n(x)"<<endl;
	cin>>n>>x;
	S=Sum(n, x);
	cout<<"Bessel function of order "<<n<<" evaluated in x= "<<x<<" is " <<S<<endl;

	
		cout<<"Finished successfully";


	getch();
	return -1;
}

[apatriarca]

E' abbastanza normale che il valore ad un certo punto converga. Il denominatore cresce molto più velocemente del numeratore e quindi il valore scende velocemente a zero. Quando la differenza tra la somma parziale e il nuovo valore sono sufficientemente grandi, il nuovo valore smetterà di avere effetto e il valore rimarrà fissato.

Supponi ad esempio di avere un numero a virgola mobile decimale fittizio con 3 cifre decimale e di avere la somma uguale a \(0.372\). Se il tuo nuovo addendo è qualcosa come \(1.73e-5,\) la somma sarà uguale a \( 0.3720173 \) che in questo sistema numerico sarà semplicemente \(0.372\) come il valore precedente della somma. Usare un numero di iterazioni uguali a 100 è in effetti in questo caso del tutto inutile.

La vera domanda da porsi è quindi un'altra: il valore a cui converge è diverso da quello che ti aspettavi? Potresti fare degli esempi? Non ho provato a compilare il codice, ma mi sembra a posto.

[irelimax]

mm si in effetti hai ragione..però stavo cercando di capire.. quando vado a calcolare la funzione di bessel per n=2 e x=3 (ad esempio) ottengo un valore 2.46604, mentre in tutte le referenze on line guardando il grafico, dovrebbe essere un valore intorno allo 0.5...inoltre quando vado a fare un grafico per le J_n in un intervallo [0,30] ad esempio, ottengo un andamento oscillante solo per n=0, mentre per n=1 e n=2 ottengo quasi delle rette...quindi ho supposto che ci fosse qualcosa che non andava e ho supposto che il problema fosse nella sommatoria..

[apatriarca]

In effetti c'è un errore. Non devi sommare \( (x/2)^n \) al resto della sommatoria, ma moltiplicarlo. Oppure puoi usarlo come valore iniziale di num e inizializzare S con num/den. Facendo il test con n=2 e x=3 ottengo 0.486091 (Wolfram Alpha restituisce lo stesso valore con lo stesso numero di cifre).

[irelimax]

Mannaggia a me! grazie mille!
Un altro dubbio (è sempre lo stesso esercizio, per questo non apro un nuovo post): devo scrivere una funzione per valutare un integrale tra 0 e infinito (di una opportuna funzione di bessel); posto che devo sostituire il valore "infinito" con un numero sufficientemente grande, pensavo di usare il metodo di simpson, ma con un passo h=0.1 (ad esempio) ho paura di appesantire il tutto.. consigliate altri metodi numerici per approssimare l'integrale o comunque come faccio a capire che il passo h=0.1 che scelgo va bene? Grazieeee

[apatriarca]

Un metodo molto semplice e che richiede relativamente poca fatica è Monte Carlo + Importance Sampling. La parte più complicata è quella di scegliere la distribuzione da usare per estrarre i campioni.

Alternativamente dovresti provare ad utilizzare un cambio di variabili in modo da avere un intervallo di integrazione finito. Puoi per esempio scrivere \( x = \tan(y)\) e integrare rispetto a \(y\) nell'intervallo \([0, \pi/2)\). Ovviamente non è l'unica trasformazione possibile..

[irelimax]

facendo un po' di ricerche sono arrivata a questa conclusione..poiche' l'integrale che mi si chiede di calcolare e' del tipo $\int_0^{+infty} e^{-ax}*J_n(bx)dx$, questo e' un caso in cui si puo' usare la formula di gausse-laguerre $\int_0^{+infty} x^ \alpha*e^{-x}*f(x)dx = \sum_0^{N-1} w_j * f(x_j)$ con $\alpha = 0$ e $f(x)$ la mia bella funzione di bessel. Adesso il mio problema e' trasferire il tutto in codice perche' non mi sembra cosi' immediato calcolare i pesi $w_j$ e le ascisse $x_j$...
suggerimenti o risorse da consulare? grazie mille
e comunque mi rimane il dubbio di quel parametro $a$ come esponente..perche' poi mi si chiede di verificare che l integrale sia uguale ad un'opportuna frazione nel caso in cui a=1 e a=2.5...quindi non so se sono nella giusta direzione..