Версия для печати темы
Нажмите сюда для просмотра этой темы в оригинальном формате
Форум программистов > C/C++: Общие вопросы > Ошибка в "Методе Симпсона"


Автор: georgjang 17.12.2008, 18:32
Код работает но в методе Симпсона есть какой то баг который мне так и не получилось выявить, перебирал его не единожды но тщетно, результат вычислений "-1.#IND". Функции в которых выявлен баг в виде коментариев. Суть кода следующая, на плоскости, фигуры- эллипс, окружность и ещё прямая их пересечением образуются сектора, площадь которых и нужно расчитать. Для наглядности задание в приложении. Просьба укажите в чём недоразумение.

(sqrt(18*18-(m+8.141)*(m+8.141)))*40/18-5 - эллипс
sqrt(24*24-(j-11.083)*(j-11.083))+2 -окружность

Код

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

void main(void)
{
    double a1,a2,a3,b1,b2,b3;
    double h1,h2,h3,e,e1;
    double j,j1,m,m1;    
    long n,i;
    int i1;
    double S[28];

    a1=0;
    a2=-12.8335;
    a3=-26;
    b1=9.68382;
    b2=0;
    b3=0;

n:    cout<<"Vvedite chislo razbieniy n = "; 
    cin >> n;
    
     h1=(b1-a1)/n;//вычисление шага интегрирования
     h2=(b2-a2)/n;
        h3=(b3-a3)/n;    
         S[4]=0; 
         S[5]=0;
         S[7]=0;
         S[6]=0;
         S[14]=0;
         S[15]=0;
         S[17]=0;
         S[16]=0;


for(i=0;i<n;i++)
{

//значени точки интегрирования

e=a1+h1*i+h1+0.5;
j=a2+h2*i+h2*.5;
m=a3+h3*i+h3*.5;

//<<<<<<<< Метод прямоугоника >>>>>>

S[1]=fabs(sqrt(24*24-(e-11.083)*(e-11.083))+2)*h1;
S[5]=S[5]+S[1]; 
S[2]=fabs(sqrt(24*24-(j-11.083)*(j-11.083))+2)*h2;
S[6]=S[6]+S[2];
S[3]=fabs((sqrt(18*18-(m+8.141)*(m+8.141)))*40/18-5)*h3;
S[7]=S[7]+S[3];  
 

//<<<<<<<< Метод трапеции >>>>>>>>

S[11]=(fabs(sqrt(24*24-((e+h1*.5)-11.083)*((e+h1*.5)-11.083))+2)+fabs(sqrt(24*24-((e-h1*.5)-11.083)*((e-h1*.5)-11.083))+2))*h1/2;
S[15]=S[15]+S[11];
S[12]=(fabs(sqrt(24*24-((j+h2*.5)-11.083)*((j+h2*.5)-11.083))+2)+fabs(sqrt(24*24-((j-h2*.5)-11.083)*(((j-h2*.5))-11.083))+2))*h2/2;
S[16]=S[16]+S[12];
S[13]=(fabs((sqrt(18*18-((m+h3*.5)+8.141)*((m+h3*.5)+8.141)))*40/18-5)+fabs((sqrt(18*18-((m-h3*.5)+8.141)*((m-h3*.5)+8.141)))*40/18-5))*h3/2;
S[17]=S[17]+S[13];

}

S[4]=S[7]-S[6];    
S[5]=S[5]-42.421;
S[4]=S[4]+S[5];

S[14]=S[17]-S[16];
S[15]=S[15]-42.421;//из общей площади сектора вычитается площадь треугольника "42.421"
S[14]=S[14]+S[15];

//<<<<<<< Метод Симпсона >>>>>><

{
        
S[25]=fabs(sqrt(24*24-(a1-11.083)*(a1-11.083))+2)+fabs(sqrt(24*24-(b1-11.083)*(b1-11.083))+2);
S[26]=fabs(sqrt(24*24-(a2-11.083)*(a2-11.083))+2)+fabs(sqrt(24*24-(b2-11.083)*(b2-11.083))+2);
S[27]=fabs((sqrt(18*18-(a3+8.141)*(a3+8.141)))*40/18-5)+fabs((sqrt(18*18-(b3+8.141)*(b3+8.141)))*40/18-5);
        
        for(i=0;i<n;i+=2)
        {

          //значени точки интегрирования

                           e1=a1+h1*i;
            j1=a2+h2*i;
            m1=a3+h3*i;

S[21]=(fabs(sqrt(24*24-((e1+h1)-11.083)*((e1+h1)-11.083))+2))*2+fabs((sqrt(24*24-((e1-h1)-11.083)*((e1-h1)-11.083))+2))*4 ;

  S[25]=S[25]+S[21];

//S[22]=(fabs(sqrt(24*24-((j1+h2)-11.083)*((j1+h2)-11.083))+2))*2+fabs((sqrt(24*24-((j1-h2)-11.083)*((j1-h2)-11.083))+2))*4 ;

    S[26]=S[26]+S[22];

//S[23]=(fabs((sqrt(18*18-((m1+h3)+8.141)*((m1+h3)+8.141)))*40/18-5))*2+fabs(((sqrt(18*18-((m1-h3)+8.141)*((m1-h3)+8.141)))*40/18-5))*4 ;

S[27]=S[27]+S[23];
S[24]=S[27]*h3/3-S[26]*h2/3;

        }

S[25]=S[25]-42.421; 
S[24]=S[24]+S[25];
}

cout<<S[4]<<'\t'<<S[14]<<'\t'<<S[24]<<endl;

cout<<"Dly prodoljeniy najmite 1 "<<endl;
cin >> i1;
     if(i1==1)
     { goto n;}
        else
         {goto end ;}
end:;
}

Автор: andrew_121 17.12.2008, 19:39
А отформатировать код слабо? Так ведь читать невозможно!

Автор: georgjang 17.12.2008, 21:31
Нет не слабо, т.к. мой статус в отличии от Вашего "новичок", думаю в праве поинтересоваться что с "этим" не так ? 

Автор: georgjang 26.12.2008, 23:12
Вот что своими стараньями смог сделать, баг остался ....... просьба "не закидывайте камнями ......."

Код

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

double Circle (double *e)

{
double S;
S=fabs(sqrt(24*24-(*e-11.083)*(*e-11.083))+2);
return (S);
}

double Circle1 (double *j)

{
double S;
S=fabs(sqrt(24*24-(*j-11.083)*(*j-11.083))+2);
return (S);
}

double Ellipse (double *m)

{
double S;
S=fabs((sqrt(18*18-(*m+8.141)*(*m+8.141)))*40/18-5);
return (S);
}

void main(void)
{
double a1,a2,a3,b1,b2,b3;
double h1,h2,h3;
double j,j1,j2,m,m1,m2,e,e1,e2; 
long n,i;
int i1;
double S[28];

a1=0;
a2=-12.8335;
a3=-26;
   b1=9.68382;
b2=0;
b3=0;

n: cout<<"Vvedite chislo razbieniy n = "; 
cin >> n;

    h1=(b1-a1)/n;//вычисление шага интегрирования
    h2=(b2-a2)/n;
 h3=(b3-a3)/n; 
     S[4]=0; 
  S[5]=0;
  S[7]=0;
  S[6]=0;
     S[14]=0;
  S[15]=0;
  S[17]=0;
  S[16]=0;

//<<<<<<<< Метод прямоугоника >>>>>>
for(i=0;i<n;i++)
{

//значени точки интегрирования

e=a1+h1*i+h1+0.5;
j=a2+h2*i+h2*.5;
m=a3+h3*i+h3*.5;

S[1]=Circle(&e)*h1;
S[5]=S[5]+S[1]; 
S[2]=Circle1(&j)*h2;
S[6]=S[6]+S[2];
S[3]=Ellipse(&m)*h3;
S[7]=S[7]+S[3];  

}

S[4]=S[7]-S[6]; 
S[5]=S[5]-42.421;
S[4]=S[4]+S[5];
e=0;
j=0;
m=0;
//<<<<<<<<<<<<>>>>>>>>>>>>>>

//<<<<<<<< Метод трапеций >>>>>>>>>

for(i=0;i<n;i++)
{
e=a1+h1*i+h1+0.5;
j=a2+h2*i+h2*.5;
m=a3+h3*i+h3*.5;

e1=e+h1*.5;
e2=e-h1*.5;
j1=j+h2*.5;
j2=j-h2*.5;
m1=m+h3*.5;
m2=m-h3*.5;

S[11]=(Circle(&e1)+Circle(&e2))*h1/2;
S[15]=S[15]+S[11];
S[12]=(Circle1(&j1)+Circle1(&j2))*h2/2;
S[16]=S[16]+S[12];
S[13]=(Ellipse(&m1)+Ellipse(&m2))*h3/2;
S[17]=S[17]+S[13];
}

S[14]=S[17]-S[16];
S[15]=S[15]-42.421;
S[14]=S[14]+S[15];
e1=0;
e2=0;
j1=0;
j2=0;
m1=0;
m2=0;
//<<<<<<<<<<<<>>>>>>>>>>>>>>

//<<<<<<< Метод Симпсона >>>>>>>

    
S[25]=Circle(&a1)+Circle(&b1);
S[26]=Circle1(&a2)+Circle1(&b2);
S[27]=Ellipse(&a3)+Ellipse(&b3);
       
for(i=0;i<n;i+=2)
{

//значени точки интегрирования

e1=a1+h1*i+h1;
j1=a2+h2*i+h2;
m1=a3+h3*i+h3;
e2=a1+h1*i-h1;
j2=a2+h2*i-h2;
m2=a3+h3*i-h3;

S[21]=(Circle(&e1))*2+(Circle(&e2))*4;    
S[25]=S[25]+S[21];
S[22]=(Circle(&j1))*2+(Circle(&j2))*4;     
S[26]=S[26]+S[22];
S[23]=(Ellipse(&m1))*2+(Ellipse(&m2))*4;
S[27]=S[27]+S[23];
S[24]=S[27]*h3/3-S[26]*h2/3;

}

S[25]=S[25]-42.421;
S[24]=S[24]+S[25];
//<<<<<<<<<<<<>>>>>>>>>>>>>>

cout<<S[4]<<'\t'<<S[14]<<'\t'<<S[24]<<endl;

cout<<"Dly prodoljeniy najmite 1 "<<endl;
cin >> i1;
    if(i1==1)
 { goto n;}
       else
  {goto end;}
end:;
}

Автор: Dmi3ev 26.12.2008, 23:36
так больше шансов разобраться... или, что кто-то станет разбираться, форматируйте всегда код...
Код

#include <math.h>
#include <iostream.h>
#include <stdio.h>
//--------------------------------------------------------------
double Circle (double *e)
{
 double S;
 S=fabs(sqrt(24*24-(*e-11.083)*(*e-11.083))+2);
 return (S);
}
//--------------------------------------------------------------
double Circle1 (double *j)
{
 double S; 
 S=fabs(sqrt(24*24-(*j-11.083)*(*j-11.083))+2);
 return (S);
}
//--------------------------------------------------------------
double Ellipse (double *m)
{
 double S; 
 S=fabs((sqrt(18*18-(*m+8.141)*(*m+8.141)))*40/18-5);
 return (S);
}
//--------------------------------------------------------------
void main(void)
{
 double a1,a2,a3,b1,b2,b3;
 double h1,h2,h3;
 double j,j1,j2,m,m1,m2,e,e1,e2;
 long n,i;
 int i1;
 double S[28];

 a1=0;
 a2=-12.8335;
 a3=-26;

 b1=9.68382;
 b2=0;
 b3=0;

 n: cout<<"Vvedite chislo razbieniy n = ";
 cin >> n;

 h1=(b1-a1)/n;
 h2=(b2-a2)/n;
 h3=(b3-a3)/n;

 S[4]=0;
 S[5]=0;
 S[7]=0;
 S[6]=0;
 S[14]=0;
 S[15]=0;
 S[17]=0;
 S[16]=0;

 for(i=0;i<n;i++)
 {
  e=a1+h1*i+h1+0.5;
  j=a2+h2*i+h2*.5;
  m=a3+h3*i+h3*.5;
  S[1]=Circle(&e)*h1;
  S[5]=S[5]+S[1];
  S[2]=Circle1(&j)*h2;
  S[6]=S[6]+S[2];
  S[3]=Ellipse(&m)*h3;
  S[7]=S[7]+S[3];
 }
 S[4]=S[7]-S[6];
 S[5]=S[5]-42.421;
 S[4]=S[4]+S[5];

 e=0;
 j=0;
 m=0;

 for(i=0;i<n;i++)
 {
  e=a1+h1*i+h1+0.5;
  j=a2+h2*i+h2*.5;
  m=a3+h3*i+h3*.5;
  e1=e+h1*.5;
  e2=e-h1*.5;
  j1=j+h2*.5;
  j2=j-h2*.5;
  m1=m+h3*.5;
  m2=m-h3*.5;
  S[11]=(Circle(&e1)+Circle(&e2))*h1/2;
  S[15]=S[15]+S[11];
  S[12]=(Circle1(&j1)+Circle1(&j2))*h2/2;
  S[16]=S[16]+S[12];
  S[13]=(Ellipse(&m1)+Ellipse(&m2))*h3/2;
  S[17]=S[17]+S[13];
 }
 S[14]=S[17]-S[16];
 S[15]=S[15]-42.421;
 S[14]=S[14]+S[15];

 e1=0;
 e2=0;
 j1=0;
 j2=0;
 m1=0;
 m2=0;


 S[25]=Circle(&a1)+Circle(&b1);
 S[26]=Circle1(&a2)+Circle1(&b2);
 S[27]=Ellipse(&a3)+Ellipse(&b3);

 for(i=0;i<n;i+=2)
  {
   e1=a1+h1*i+h1;
   j1=a2+h2*i+h2;
   m1=a3+h3*i+h3;
   e2=a1+h1*i-h1;
   j2=a2+h2*i-h2;
   m2=a3+h3*i-h3;
   S[21]=(Circle(&e1))*2+(Circle(&e2))*4;
   S[25]+=S[21];
   S[22]=(Circle(&j1))*2+(Circle(&j2))*4;
   S[26]+=S[22];
   S[23]=(Ellipse(&m1))*2+(Ellipse(&m2))*4;
   S[27]+=S[23];
   S[24]=S[27]*h3/3-S[26]*h2/3;
  }
 S[25]+=-42.421;
 S[24]+=S[25];
 //<<<<<<<<<<<<>>>>>>>>>>>>>>
 cout<<S[4]<<'\t'<<S[14]<<'\t'<<S[24]<<endl;
 cout<<"Dly prodoljeniy najmite 1 "<<endl;
 cin >> i1;
 if(i1==1)
  goto n;
 else
  goto end;
end:;
}

а вот это можно и покороче записывать...
Цитата

double Circle (double *e)
{
 double S;
 S=fabs(sqrt(24*24-(*e-11.083)*(*e-11.083))+2);
 return (S);
}

вот так:
Код

double Circle (double *e) {
 return (fabs(sqrt(24*24-(*e-11.083)*(*e-11.083))+2));
}

Powered by Invision Power Board (http://www.invisionboard.com)
© Invision Power Services (http://www.invisionpower.com)