【實驗性質】綜合性實驗。

【實驗目的】理解插值型積分法；掌握復化積分法算法。?

【實驗內容】

1對 ，用復化梯形積分和變步長梯形積分求值（截斷誤差不超過）。

【理論基礎】

積分在工程中有重要的應用，數值積分的基本思想是用被積函在區間上的一些點處的值 的線性組合作為積分的近似值：

實際應用中 f x( )是未知的，一般用 f x( ) 的次數不超過n的插值多項式來代替（插值型求積方法）：

【實驗過程】

1.用復化梯形積分求解，給出代碼，并用表格記載求解過程〔區間數n=30、50、70、100、150〕。

程序代碼：

頭文件：

#ifndef DEFINITEINTEGRAL_H

#define DEFINITEINTEGRAL_H

#include<math.h>

class definiteintegral

{

public:

??? definiteintegral();

??? double T_Integral(double a ,double b,double(*f)(double));

??? double S_Integral(double a ,double b,double(*f)(double));

??? double C_Integral(double a ,double b,double(*f)(double));

??? double Tn_Integral(double a ,double b,double(*f)(double),int n);

??? double Sn_Integral(double a ,double b,double(*f)(double),int n);

??? double Cn_Integral(double a ,double b,double(*f)(double),int n);

??? double V_T_Integral(double a, double b, double(*f)(double),double e);

};

#endif // DEFINITEINTEGRAL_H

主函數：

//實驗六

#include <iostream>

#include <windows.h>

#include "colvector.h"

#include "matrix.h"

#include <windows.h>

#include "linearequations.h"

#include "interpolationpolynomial.h"

#include "definiteintegral.h"

double f1(double x){

??? return x*x*x-sin(x)-4*x-1;

}

double f2(double x){

??? return 3*x*x-cos(x)-4;

}

double f3(double x){

??? double result =sin(x)+4*x-1;

??????? if(result>0){

??????????? return pow(result,1.0/3);

??????? }else{

???????????? return -pow(fabs(result),1.0/3);

??????? }

}

using namespace std;

double f4(double x)

{

??? return 1/(1+x*x);

}

int main()

{

??? SetConsoleOutputCP(CP_UTF8);

??? double a=0,b=1.0;

??? definiteintegral obj;

??? cout<<"梯形積分:\t\t"<<obj.T_Integral(a,b,f4)<<endl;

??? cout<<"Simpson積分:\t\t"<<obj.S_Integral(a,b,f4)<<endl;

??? cout<<"Cotess積分:\t\t"<<obj.C_Integral(a,b,f4)<<endl;

??? cout<<endl;

??? cout<<"復化梯形積分:\t\t"<<obj.Tn_Integral(a,b,f4,150)<<endl;

??? cout<<"復化Simpson積分:\t"<<obj.Sn_Integral(a,b,f4,30)<<endl;

??? cout<<"復化Cotes積分:\t\t"<<obj.Cn_Integral(a,b,f4,30)<<endl;

??? cout<<endl;

??? cout<<"交步長梯形積分:\t\t"<<obj.V_T_Integral(a,b,f4,0.00001)<<endl;

??? return 0;

}

代碼塊：

#include "definiteintegral.h"

definiteintegral::definiteintegral(){}

double definiteintegral::T_Integral(double a,double b,double(*f)(double)){

??? return (b-a)/2.0*(f(a)+f(b));

}

double definiteintegral::S_Integral(double a,double b,double(*f)(double)){

??? return (b-a)/6.0*(f(a)+4*f((a+b)/2)+f(b));

}

double definiteintegral::C_Integral(double a,double b,double(*f)(double)){

??? double h=(b-a)/4.0;

??? return (b-a)/90*(7*f(a)+32*f(a+h)+12*f(a+2*h)+32*f(a+3*h)+7*f(b));

}

double definiteintegral::Tn_Integral(double a,double b,double(*f)(double),int n){

??? double h=(b-a)/n;

??? double result =0;

??? double sum=0;

??? for(int i=1;i<=n-1;i++){

??????????? sum +=f(a+i*h);

??? }

??? result =h/2.0*(f(a)+2*sum+f(b));

??? return result;

}

double definiteintegral::Sn_Integral(double a,double b,double(*f)(double),int n){

??? double h=(b-a)/n;

??? double result =0;

??? double sum1=0,sum2=0;

??? for(int i=1;i<=n;i++){

??????????? sum1 +=f(a+(i-1)*h+h/2);

??? }

??? for(int i=1;i<=n-1;i++){

??????????? sum2 +=f(a+i*h);

??? }

??? result =h/6.0*(f(a)+4*sum1+2*sum2+f(b));

??? return result;

}

double definiteintegral::Cn_Integral(double a,double b,double(*f)(double),int n){

??? double h=(b-a)/n;

??? double result =0;

??? double sum1=0,sum2=0,sum3=0,sum4=0;

??? for(int i=1;i<=n;i++){

??????????? sum1 +=f(a+(i-1)*h+h/4);

??? }

??? for(int i=1;i<=n;i++){

??????????? sum2 +=f(a+(i-1)*h+h/2);

??? }

??? for(int i=1;i<=n;i++){

??????????? sum3 +=f(a+(i-1)*h+3*h/4);

??? }

??? for(int i=1;i<=n-1;i++){

??????????? sum4 +=f(a+i*h);

??? }

??? result =h/90.0*(7*f(a)+32*sum1+12*sum2+32*sum3+14*sum4+7*f(b));

??? return result;

}

double definiteintegral::V_T_Integral(double a, double b, double(*f)(double),double e){

??? double h=(b-a);

??? double T1=h/2.0*(f(a)+f(b));

??? double T2=T1/2.0+h/2.0*f(a+h/2.0);

??? double error=fabs(T2-T1);

??? while(error>=e){

??????? T1=T2;

??????? h=h/2;

??????? double x=a+h/2;

??????? double sum=0;

??????? while(x<b){

??????????? sum+=f(x);

??????????? x+=h;

??????? }

??????? T2=T1/2.0+h/2.0*sum;

??????? error=fabs(T2-T1);

??? }

??? return T2;

}

表格：

n=30

n=50

n=70

n=100

n=150

0.785352

0.785381

0.785390

0.785394

0.785396

?

2.用變步長梯形積分求解，給出代碼，并用表格記載求解過程。

代碼塊：

double definiteintegral::V_T_Integral(double a, double b, double(*f)(double),double e){

??? double h=(b-a);

??? double T1=h/2.0*(f(a)+f(b));

??? double T2=T1/2.0+h/2.0*f(a+h/2.0);

??? double error=fabs(T2-T1);

??? while(error>=e){

??????? T1=T2;

??????? h=h/2;

??????? double x=a+h/2;

??????? double sum=0;

??????? while(x<b){

??????????? sum+=f(x);

??????????? x+=h;

??????? }

??????? T2=T1/2.0+h/2.0*sum;

??????? error=fabs(T2-T1);

??? }

??? return T2;

}

表格：

3.比較復化積分與變步長積分，分析實驗出現的問題，總結解決辦法。

??? 復化積分是將一個區間分成若干子區間，然后在每個子區間上應用數值積分方法。它的優點是簡單易實現，計算結果比較穩定。但是，如果子區間的數量不夠多，或者函數在某些子區間上變化較大，可能會導致計算結果的誤差較大。為了解決這個問題，可以增加子區間的數量，或者使用自適應方法，根據函數的變化情況來調整子區間的數量。

變步長積分是根據函數的變化情況，調整積分步長來提高計算精度。它的優點是能夠更好地適應函數的變化情況，減小誤差。但是，如果調整步長的策略不合理，可能會導致計算時間過長。為了解決這個問題，可以采用適當的步長調整策略，如自適應選取步長或者根據函數的一階或二階導數來調整步長。

【實驗心得】

在本次實驗中，我們學習了三種數值積分的方法，包括差值積分法、復化積分法和復化梯形公式。通過這些實驗，我對數值積分的原理和計算方法有了更深入的了解。

在差值積分法中，我們使用了牛頓-科特斯公式對函數進行了差值，然后通過對差值多項式進行求和來計算積分。這種方法的優點是計算簡單，適用于低次多項式的積分，但對于高次多項式的積分誤差較大。

復化積分法是將計算區間分成若干小區間，然后對每個小區間應用數值積分方法。我們采用了復化梯形公式，其原理是通過將每個小區間近似為梯形來計算積分。這種方法誤差較差值積分法要小，但計算較復雜。

綜合考慮精度和計算復雜度，我們可以選擇合適的數值積分方法。如果函數是低次多項式，可以使用差值積分法進行計算。對于復雜函數，可以采用復化積分法進行分區間計算。通過這些實驗，我掌握了數值積分的基本原理和計算方法，并且了解了不同方法的優缺點，這對于解決實際問題具有重要的參考價值。

得??? 分_____________

?

評閱日期_____________

?

教師簽名_____________