Метод Ньютона Котеса

[samuel21]

Для решения определенных интегралов необходимо написать метод Ньютона Котеса.
Вот что уже готово:

PHP код:

      double f_Ex(double x){return x*x/9;}

      void newton(){
      int i,j,n=10,N=5;
      double b=5,a=0,S,h,A[N],X[n],x[N],Y[N],Sum,dh,**H;
 
       H=new double*[5];
       for(i=0;i<6;i++) H[i]=new double[6];




     h=(b-a)/n;
     X[0]=a;
     
     for(i=1;i<n;i++)X[i]=X[i-1]+h; 
     
     Sum=0; 
     for(i=0;i<n-1;i++){ 
                        dh=(X[i+1]-X[i])/N;
                        x[0]=X[i];
                        for(j=1;j<N;j++)x[j]=x[0]+j*dh;
                        for(j=0;j<N;j++)Y[j]=f_Ex(x[j]);
                        S=0;
                        for(j=0;j<N;j++)A[j]=(X[i+1]-X[i])*H[N][j];
                        for(j=0;j<N;j++)S+=Y[j]*A[j];
                        Sum+=S;
                        }
    printf("%lf\n",Sum);                          
     
} 


Результат ~8,9, а должно быть ~4,8. Что не правильно?
Делаю по данной блок схеме:

Содержимое матрицы H (постоянные Котеса)

PHP код:

H[0][0]=1/2;
H[0][1]=1/2;

H[1][0]=1/3;
H[1][1]=4/3;
H[1][2]=1/3;

H[2][0]=3/8;
H[2][1]=9/8;
H[2][2]=9/8;
H[2][3]=3/8;

H[3][0]=14/45;
H[3][1]=64/45;
H[3][2]=8/15;
H[3][3]=64/45;
H[3][4]=14/45;

H[4][0]=95/288;
H[4][1]=125/96;
H[4][2]=125/144;
H[4][3]=125/144;
H[4][4]=125/96;
H[4][5]=95/288;

H[5][0]=41/140;
H[5][1]=54/35;
H[5][2]=27/140;
H[5][3]=68/35;
H[5][4]=27/140;
H[5][5]=54/35;
H[5][6]=41/140; 


И вопрос по самому методу: есть ли формулы вычисления постоянных (не в интегральном виде), что бы не заполнять матрицу вручную?

[Ezhik Kihze]

H=new double*[5]; 
       for(i=0;i<6;i++) H[i]=new double[6];

Нету у H по первому индексу 5, а по второму 6... Выход за пределы массива.

     for(i=0;i!=n;i++)X[i]=X[i-1]+h; 
      
     Sum=0; 
     for(i=0;i!=n-1;i++){

У вас размерность массивов 5, а бегаете как по массиву с размерностью 10

[samuel21]

У вас размерность массивов 5, а бегаете как по массиву с размерностью 10

По массиву H(5) я бегаю тут:

PHP код:

for(j=0;j<N;j++)A[j]=(X[i+1]-X[i])*H[N][j]; 


Первый недочет исправил, но разницы нет. Не работает

[samuel21]

Ошибки были в циклах, долго возился, результат:

PHP код:

h=(b-a)/n;
     X[0]=a; 
     for(i=1;i<n;i++)X[i]=X[i-1]+h;
     Sum=0; 

     for(i=0;i<n;i++){ 
                        if(i==0)dh=(X[i+1]-X[i])/N;
                        else dh=(X[i]-X[i-1])/N;
                        x[0]=X[i]; 
                        for(j=0;j<N+1;j++)x[j]=x[0]+(j*dh);
                        for(j=0;j<N+1;j++)y[j]=f_Ex(x[j]);
                        S=0; 
                        for(j=0;j<N+1;j++){
                                           if(i==0) A[j]=(X[i+1]-X[i])*H[N][j];
                                           else A[j]=(X[i]-X[i-1])*H[N][j];
                                           }
                        for(j=0;j<N+1;j++)S=S+y[j]*A[j];
                        Sum+=S;
                      }