Числовой интеграл от 0 до бесконечности

Question

Числовой интеграл от 0 до бесконечности

Моя цель состоит в том, чтобы вычислить числовой интеграл функции распределения вероятности (PDF) расстояния электрона от ядра атома водорода на языке программирования Си. Я написал пример кода, однако он не может правильно найти числовое значение из-за того, что я не могу увеличить предел настолько, насколько это необходимо, по моему мнению. Я также включил библиотеку, но не могу использовать значения, указанные в следующем посте, в качестве целых границ: минимальное и максимальное значение типа данных в Си. Какое лекарство в этом случае? Может стоит переключиться на другой язык программирования? Любая помощь и предложение приветствуется, спасибо заранее.

Изменить: После некоторого значения я получаю ошибку сегментации ошибки. Я проверил фактический результат интеграла, чтобы быть 0.0372193 с Wolframalpha. В дополнение к этому, если я увеличиваю k в меньших количествах, в результате я получаю ноль, поэтому я определил r[k]=k, я знаю, что он должен быть меньше для повышения точности.

#include <stdio.h>
#include <math.h>
#include <limits.h>
#define a0 0.53 
int N = 200000;
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[], long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {

    long double P[N], r[N], a;
    // Declare and initialize the loop variable
    int k = 0;
    for (k = 0; k <  N; k++)
    {
        r[k] = k ;
        P[k] = r[k] * r[k] * exp( -2*r[k] / a0);
        //printf("%.20Lf \n", r[k]);
        //printf("%.20Lf \n", P[k]);
    }
    a = trapezoid(r, P);    
    printf("%.20Lf \n", a);
}

Последний код:

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#define a0 0.53 
#define N LLONG_MAX
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[],long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {
    printf("%Ld", LLONG_MAX);
    long double * P = malloc(N * sizeof(long double));
    long double * r = malloc(N * sizeof(long double));
    // Declare and initialize the loop variable
    int k = 0;
    long double integral;
    for (k = 1; k <  N; k++)
        {
        P[k] = r[k] * r[k] * expl( -2*r[k] / a0);
        }
    integral = trapezoid(r, P);
    printf("%Lf", integral);

}

Редактировать последний код работает:

#include <stdio.h>
#include <math.h>
#include <limits.h>
#include <stdlib.h>
#define a0 0.53 
#define N LONG_MAX/100
// This value of N is the highest possible number in long double
// data format. Change its value  to adjust the precision of integration
// and computation time.
// The discrete integral may be defined as follows:
long double trapezoid(long double x[],long double f[]) {
  int i;
  long double dx = x[1]-x[0];
  long double sum = 0.5*(f[0]+f[N]);

    for (i = 1; i <  N; i++) 
        sum+=f[i];
  return sum*dx;  
}
main() {
    printf("%Ld \n", LLONG_MAX);
    long double * P = malloc(N * sizeof(long double));
    long double * r = malloc(N * sizeof(long double));
    // Declare and initialize the loop variable
    int k = 0;
    long double integral;
    for (k = 1; k <  N; k++)
        {
        r[k] = k / 100000.0;
        P[k] = r[k] * r[k] * expl( -2*r[k] / a0);
        }
    integral = trapezoid(r, P);
    printf("%.15Lf \n", integral);
    free((void *)P);
    free((void *)r);
}

В частности, я изменил определение для r [k], используя число с плавающей точкой в операции деления для получения длинного двойного в результате, а также, как я уже говорил в своем последнем комментарии, я не могу перейти на Ns больше, чем LONG_MAX/100 и я думаю, что я должен исследовать код и malloc дальше, чтобы получить проблему. Я нашел точное значение, которое получается аналитически, принимая пределы; Я подтвердил результат с TI-89 Titanium и Wolframalpha (как в численном, так и в аналитическом отношении), не считая того, чтобы делать это самому. Правило трапеции сработало довольно хорошо, когда размер интервала был уменьшен. Большое спасибо за все плакаты здесь за их идеи. Имея значение 2147483647, LONG_MAX не так уж и велик, как я ожидал, кстати, не должен ли предел быть около десяти до 308?

1

c gcc ubuntu-12.04 numerical-integration

Источник

user2840680 01 ноя '13 в 09:58

1 ответ

Решение

Другие вопросы по тегам c gcc ubuntu-12.04 numerical-integration

user1139697 01 ноя '13 в 10:11 2013-11-01 10:11 · Accepted Answer · 2013-11-01 10:11

Численная точка зрения

Обычный метод трапеции не работает с неправильными интегралами. Таким образом, квадратурные правила Гаусса намного лучше, поскольку они не только обеспечивают точность 2n-1 (то есть для многочлена степени 2n-1 они возвращают правильное решение), но и управляют неправильными интегралами, используя правильную весовую функцию.,

Если ваш интеграл неправильный с обеих сторон, вы должны попробовать квадратуру Гаусса-Эрмита, в противном случае используйте квадратуру Гаусса-Лагерра.

Ошибка "переполнения"

long double P[N], r[N], a;

P имеет размер примерно 3 МБ, и поэтому r, Это слишком много памяти. Выделите память вместо:

long double * P = malloc(N * sizeof(long double));
long double * r = malloc(N * sizeof(long double));

Не забудьте включить <stdlib.h> и использовать free на обоих P а также r если они тебе больше не нужны. Кроме того, вы не можете получить доступ к N-й записи, поэтому f[N] неправильно.

Использование квадратуры Гаусса-Лагерра

Теперь Гаусс-Лагерр использует exp(-x) как весовая функция. Если вы не знакомы с квадратурой Гаусса: результат E(f) является интегралом w * f, где w это весовая функция.

Ваш f выглядит так, и:

f x = x^2 * exp (-2 * x / a)

Подождите минуту. f уже содержит exp(-term)так что мы можем заменить х t = x * a /2 и получить

f' x = (t * a/2)^2 * exp(-t) * a/2

поскольку exp(-t) это уже часть нашей весовой функции, теперь ваша функция идеально вписывается в квадратуру Гаусса-Лагерра. Полученный код

#include <stdio.h>
#include <math.h>

/* x[] and a[] taken from 
 * https://de.wikipedia.org/wiki/Gau%C3%9F-Quadratur#Gau.C3.9F-Laguerre-Integration
 * Calculating them by hand is a little bit cumbersome
*/
const int gauss_rule_length = 3;
const double gauss_x[] = {0.415774556783, 2.29428036028, 6.28994508294};
const double gauss_a[] = {0.711093009929, 0.278517733569, 0.0103892565016};

double f(double x){
    return x *.53/2 * x *.53/2 * .53/2;
}

int main(){
    int i;
    double sum = 0;
    for(i = 0; i < gauss_rule_length; ++i){
        sum += gauss_a[i] * f(gauss_x[i]);
    }
    printf("%.10lf\n",sum); /* 0.0372192500 */
    return 0;
}