龙贝格算法.cpp

来自「数值分析（计算方法）：欧拉方法、高斯－赛德尔迭代法、拉格朗日及分段线性插值、龙贝」· C++ 代码 · 共 88 行

#include<iostream>
#include<cmath>
#include <iomanip>
using namespace std;

double f(double x)
{  if(x==0) return 1;
   else return sin(x)/x;
}
    
int main()
{  int i=0;
   double r[20];
   double a,b,h,x,k,t1,t2,s,s1,s2,c1,c2,r1,r2;
   long double e;   
   cout<<"input a:";
   cin>>a;
   cout<<"input b:";
   cin>>b;
   cout<<"input e:";
   cin>>e;
   h=b-a;
   t1=h*(f(a)+f(b))/2;
   r[i]=t1;
   i++;                     
   k=1;
for(;;)
	{  s=0;x=a+h/2;
     do{ s=s+f(x);
         x=h+x;
       }while(x<b);
     t2=t1/2+(h/2)*s;
     r[i]=t2;i++;                
     s2=t2+(t2-t1)/3;
     r[i]=s2;i++;               
     
     if(k==1)
     {   k=k+1;h=h/2;
         t1=t2;s1=s2;
         continue;
     }
     
     c2=s2+(s2-s1)/15.0;
     r[i]=c2;i++;                   
     
     if(k==2)
     {   c1=c2;
         k=k+1;h=h/2;
         t1=t2;s1=s2;
         continue;
     }
     
     r2=c2+(c2-c1)/63.0;
     r[i++]=r2;
     
     if(k==3)
     {   r1=r2;
         c1=c2;
         k=k+1;h=h/2;
         t1=t2;s1=s2;
         continue;
     } 
     
     r1=r2;
     c1=c2;
     k=k+1;h=h/2;
     t1=t2;s1=s2;
     
     if(abs(r2-r1)<e)
        break;
        
  }
  
   int l=12,n=13;
   cout<<" k "<<setw(6)<<"T"<<setw(l)<<"S"<<setw(l)<<"C"<<setw(l)<<"R"<<endl;
   cout<<" 0 "<<setw(4)<<showpoint<<setprecision(7)<<r[0]<<endl;
   cout<<" 1 "<<setw(4)<<r[1]<<setw(n)<<r[2]<<endl;
   cout<<" 2 "<<setw(4)<<r[3]<<setw(n)<<r[4]<<setw(n)<<r[5]<<endl;
   cout<<" 3 "<<setw(4)<<r[6]<<setw(n)<<r[7]<<setw(n)<<r[8]<<setw(n)<<r[9]<<endl;
   cout<<" 4 "<<setw(4)<<r[10]<<setw(n)<<r[11]<<setw(n)<<r[12]<<setw(n)<<r[13]<<endl;
   return 0;
}