Численное интегрирование с правилом трапеции между 0 и бесконечностью Sin (x) / sqrt (x) с программированием C - PullRequest
Новая
       18

Численное интегрирование с правилом трапеции между 0 и бесконечностью Sin (x) / sqrt (x) с программированием C

1 голос
/ 18 марта 2020

У меня проблема с поиском хорошего приближения к интегралу между 0 и бесконечностью sin (x) / sqrt (x) в C программировании. Я пытаюсь использовать правило трапеции. В моем коде я также хочу, чтобы пользователь ввел значение точности, такое как 0,001, где выходное значение будет с точностью до этого количества десятичных разрядов. Что не так с этим кодом? 

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

int main(void)
{
int i, N; // integers to interate on in loops //
double h, x, y, precision, lowerLim = 0.0001, upperLim = 10000, f0, fN;

printf("accuracy:");
scanf("%lf", &precision);
N = 10;
double areastored, newarea; 
do {
    newarea = 0.0;
    areastored = newarea; // areastored is the area that I want to compare to the new area calulated as the N (number of partitions) to check the precision of the new value to see if it a better approximation //
    h = (upperLim - lowerLim)/(N-1);
    fN = sin(upperLim)/sqrt(upperLim);
    f0 = sin(lowerLim)/sqrt(lowerLim); // end points evaluated in function //
    newarea = newarea + 0.5*h*(f0 + fN); 

    for (int i = 1; i < N; i++) {
    x = lowerLim + h*i;
    y = sin(x)/sqrt(x);
    newarea = newarea + y*h;   // this loop adds all the middle trapezia areas //
    }

    printf("at N %d integral %f\n", N, newarea);
    N = N*5; // iterate the N so next round of the loop it will approximate an area with more and smaller trapezia //
} while ( fabs ( newarea - areastored ) > precision ); // if this is false then should have an area to the desired precision //

    printf("The integral evaluates to: %lf\n", newarea); 
}

Проблема в том, что, если я введу точность 0,01, область рассчитывается для N = 10, 50, 250, но не может продолжаться, а последняя область = 8.53, что не соответствует значению 1.253 ... Я ожидаю

РЕДАКТИРОВАТЬ: я внес предлагаемые изменения в приведенный выше код, что можно увидеть сейчас благодаря комментариям нескольких пользователей ниже. Спасибо большое за помощь! Теперь у меня проблема с выводом, см. Прикрепленное изображение моего терминала. Он должен был остановиться на 9-й итерации N для этого ввода, и поэтому напечатанная вейл довольно озадачивает. Почему это произошло? Еще раз спасибо за любую помощь заранее!

enter image description here

Ответы [ 2 ]

1 голос
/ 18 марта 2020

В вашем коде есть как минимум две ошибки 

// ...
newarea = 0.0;                       // <-- It's initialized here
do {
    // ...
    h = (upperLim - lowerLim)/N-1;   // <-- This doesn't do what you think
    // ...
    // It's updated, but never reset 
    newarea = newarea + 0.5*h*(f0 + fN);
    // ...
} while ( ... );

Вы должны переместить строку newarea = 0.0; внутри l oop и изменить формулу, которая вычисляет h.

Также обратите внимание, что у вас есть фиксированный верхний предел (1000), в то время как для этого типа интегралов вы должны учитывать увеличивающийся верхний предел и, возможно, неравномерный интервал сетки.

0 голосов
/ 18 марта 2020

Если вы отладите свою программу, вы увидите, что 

h = (upperLim - lowerLim)/N-1;

приводит к отрицательному значению h для N=1250.

Это приводит к бесконечному l oop 

for (x=lowerLim+h; x < upperLim; x+=h) {

, поскольку, когда x достигает достаточно большого абсолютного значения, оно больше не изменяется при добавлении h.

Вы, вероятно, имеете в виду 

h = (upperLim - lowerLim) / (N-1);
...