Моя цель — рассчитать численный интеграл функции распределения вероятностей (PDF) расстояния электрона от ядра атома водорода на языке программирования C. Я написал пример кода, однако он не может правильно найти числовое значение из-за того, что я не могу увеличить предел настолько, насколько это необходимо, по моему мнению. Я также включил библиотеку, но не могу использовать значения, указанные в следующем сообщении, в качестве интегральных границ: минимальное и максимальное значение типа данных в C . Какое средство в этом случае? Может стоит перейти на другой язык программирования? Любая помощь и предложение приветствуются, спасибо заранее.
Изменить: после некоторого значения я получаю ошибку сегментации ошибки. Я проверил, что фактический результат интеграла равен 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?
