數(shù)值分析：數(shù)值積分

一、前言

本文主要討論定積分的數(shù)值解箕般，即數(shù)值積分耐薯。為什么需要數(shù)值積分：假設 $f(x)$ 在區(qū)間上可積， $F(x)$ 為其原函數(shù)丝里；現(xiàn)實中的多數(shù)情況是：

原函數(shù) $F(x)$ 無法用初等函數(shù)表示曲初，即找不到原函數(shù)；
原函數(shù) $F(x)$ 可以找到杯聚，但是求解過程太復雜太困難了臼婆。

上面兩種情況其實非常常見，因此我們就需要用數(shù)值積分幌绍！找原定積分的"近似解"颁褂；只要精度夠高，近似解完全可以替代真實解傀广。

數(shù)值積分的方法可以分為兩大類：拉格朗日型(單步法)痢虹、高斯型(多步法)。

拉格朗日型：梯形求積主儡、辛普森求積奖唯、牛頓-科茨求積、復化求積糜值、龍貝格求積丰捷；
高斯型：高斯-勒讓德求積坯墨、高斯-拉蓋爾求積、高斯-埃爾米特求積病往。

所有數(shù)值積分的核心思想：用一個積分易求的多項式去近似替代被積函數(shù)捣染，然后對這個多項式再求定積分。 這種思想是基于插值理論停巷。

本文所有提到的方法耍攘，完整matlab程序出處在這里。

二畔勤、拉格朗日型數(shù)值積分

拉格朗日型數(shù)值積分蕾各，是采用各階的拉格朗日插值多項式來近似被積函數(shù) $f(x)$ 的。所以庆揪，這種方法要想提高精度式曲，只能增多插值節(jié)點(或者將積分區(qū)間劃分的份數(shù)增多)。和插值理論一樣缸榛，插值節(jié)點可以是隨意的吝羞！一般我們?yōu)榱朔奖憔幊滩胚x擇插值節(jié)點是各個等分點。

拉格朗日型數(shù)值積分可以再細分為兩大類：1. 基礎方法内颗；2. 進階方法钧排。基礎方法：有多少插值節(jié)點均澳，就用幾階多項式去插值近似恨溜；進階方法：分段線性插值，即不管插值節(jié)點再多负懦，相鄰2點間用線性拉格朗日插值代替筒捺！總體用分段函數(shù)近似即可柏腻。

基礎方法：梯形求積公式纸厉、辛普森求積公式、牛頓-科茨求積公式五嫂。
進階方法：復化梯形公式颗品、復化辛普森公式、復化梯形公式加密沃缘、龍貝格公式躯枢。

2.1 基礎方法

（1）梯形求積公式：積分區(qū)間 $[a,b]$ ，2個插值節(jié)點為(頭和尾)：

$x_0 = a槐臀，x_1 = b锄蹂；y_0 = f(x_0)，y_1 = f(x_1)$

用2點一階線性拉格朗日插值多項式代替被積函數(shù)并帶入原式可得：

$\int_{a}^水慨f(x)dx ≈ \frac{h}{2}[f(x_0) + f(x_1)] = \frac{b-a}{2}[f(a) + f(b)]$

（2）辛普森求積公式：將積分區(qū)間 $[a,b]$ 二等分得糜，取 $h = \frac{b-a}{2}$ 敬扛；3個插值節(jié)點為：

$x_0 = a，x_1 = \frac{b+a}{2}朝抖，x_2 = b啥箭；y_0 = f(x_0)，y_1 = f(x_1)治宣，y_2 = f(x_2)$

用3點二次拉格朗日插值多項式代替被積函數(shù)并帶入原式可得：

$\int_{a}^急侥f(x)dx ≈ \frac{b-a}{6}[f(a) + 4f\frac{a+b}{2} + f(b)]$

（3）牛頓-科茨求積公式：將積分區(qū)間 $[a,b]$ n等分，取 $h = \frac{b-a}{n}$ 侮邀；共有n+1個插值節(jié)點坏怪。注意一點：第一個節(jié)點是 $x_0$ ，是區(qū)間左端點豌拙；最后一個節(jié)點是 $x_n$ 陕悬，是區(qū)間的右端點。

$\int_{a}^按傅f(x)dx ≈ (b-a)\sum_{i=0}^{n}C_{i}^{(n)}f(x_i)$

$C_{i}^{(n)} = \frac{(-1)^{n-i}}{n\cdot i!(n-i)! }\int_{0}^{n}(t-0)(t-1)\cdots[t-(i-1)][t-(i+1)][t - (i+1)]\cdots (t-n)dt$

公式說明：i是節(jié)點的下標捉超，范圍是 $[0,n]$ ，上式連乘中i的變化是從0增加到n唯绍！一共n次循環(huán)拼岳，每次i都不能等于一個值！上式中求 $C_{i}^{(n)}$ 其實是一個很長的公式况芒！用語言不好形容惜纸，直接看下面的程序：

% n點牛頓-科茨求積公式
% 思路: 先將[a,b]區(qū)間n等分, 然后n+1個點用n+1點(n次)拉格朗日插值
% 缺點: 高階拉格朗日插值對于復雜彎曲函數(shù)龍貝格現(xiàn)象嚴重！不要插過6階绝骚！
% 以f(x) = 4/(1+x^2)函數(shù)為例

clear; clc;

up = double(input('輸入積分上限:'));
low = double(input('輸入積分下限:'));
n = double(input('將積分區(qū)域幾等分(幾階拉格朗日插值多項式):'));

% 算當前n等分牛頓-科茨系數(shù)Cn 
Cn = zeros(1,n+1); 
syms t;
% s=sym(1)符號變量類型的1, 主要是為了一開始的循環(huán)相乘
for i = 0:n
    s = sym(1);
    for num = 0:n
        if num ~= i
            s = s*(t-num);
        end
    end
    Cn(i+1) = (-1)^(n-i)/( n*factorial(i)*factorial(n-i) )*int(s,t,0,n)
end

syms x;
f = 4/(1+x^2);  % 每次修改這里的函數(shù)
R = int(f,x,low,up);    % 真實結(jié)果
x = low:(up-low)/n:up;  % n等分后每個點的x值
fx = double(subs(f));   % n等分后每個x點對應的函數(shù)值

% 近似結(jié)果為:
NC = (up-low)*sum(Cn.*fx);
fprintf('牛頓-科茨求積公式近似結(jié)果為: %.7f\n',NC);
fprintf('真實結(jié)果為: %f\n',R);

2.2 進階方法

（1）復化梯形公式：區(qū)間劃分和牛頓-科茨一樣耐版，只是復化梯形采用的是分段線性插值，即相鄰兩插值節(jié)點用線性函數(shù)代替压汪，總體用分段函數(shù)近似即可粪牲。插值節(jié)點為：

$x_i = a + ih \quad (i = 0,1,2,\cdots, n)$

n是區(qū)間等分數(shù)靴寂，h為每一段長葱轩。近似計算公式為：

$\int_{a}^f(x)dx ≈ \frac{b-a}{2n}[f(a) + 2\sum_{i=1}^{\color{red}{n-1}}f(a + i\cdot\frac{b-a}{n}) + f(b)]$

復化梯形公式是精度優(yōu)秀方法中最簡單于置、最萬能穿香、最推薦使用的方法亭引。

% n點復化求積公式——萬能！
% 思路: 先將[a,b]區(qū)間n等分, 然后每兩點間用線性拉格朗日插值——分段線性插值而已皮获。
% 優(yōu)勢: 避免高階牛頓-科茨(高階)出現(xiàn)龍貝格現(xiàn)象, 因此等分越多越好!
% 以f(x) = 4/(1+x^2)函數(shù)為例

clear; clc;

syms x;
f = 4/(1+x^2);   % 每次這里修改不同函數(shù)即可

up = double(input('輸入積分上限:'));
low = double(input('輸入積分下限:'));
n = double(input('將積分區(qū)間幾等分:'));
% 真實定積分的結(jié)果: 
R = int(f,x,low,up); 

xnum = low:(up-low)/n:up;
fmiddle = zeros(1,n-1);  % 對中間那些結(jié)果的記錄
for num = 1:n-1
    x = low + num*(up-low)/n;
    fmiddle(num) = double(subs(f));
end

x = up;
fup = double(subs(f));
x = low;
flow = double(subs(f));
% 復化梯形公式近似結(jié)果:
FT = (up-low)/(2*n)*( flow + 2*sum(fmiddle) + fup ); 

fprintf('復化梯形求積公式近似結(jié)果為: %.7f\n',FT);
fprintf('真實結(jié)果為: %f\n',R);

（2）復化辛普森公式：用兩種理解方式：

將區(qū)間分成偶數(shù)份焙蚓，非重疊的相鄰3點間用一個二次拉格朗日插值代替；
將區(qū)間分成n等分，再在每一份中插入一個中心節(jié)點购公；這樣每一份其實里面就有3個點赵哲。

無論用哪種方式，都是分段2次多項式插值而已君丁。本文推薦采用第一種方式：將區(qū)間[a,b]等分為偶數(shù)份n = 2m枫夺，每個小區(qū)間長度h = (b-a)/(2m)，插值節(jié)點為：

$x_i = a + ih \quad (i = 0,1,2,\cdots, 2m)$

每3個非重疊相鄰點 $[x_{2i-1},x_{2i}]$ 用二次拉格朗日插值代替即可的最后公式為：

$\int_{a}^绘闷f(x)dx ≈ \frac{b-a}{6m}[f(a) + 4\sum_{i=1}^{\color{red}{m}}f(a + (2i-1)\cdot\frac{b-a}{2m}) + f(b)]$

公式若不太清楚關(guān)系可直接看程序：

% n點復化辛普森公式——萬能橡庞！比復化梯形更準！
% 思路: 先將[a,b]區(qū)間n等分, 然后每兩點間再插入一個中心節(jié)點印蔗！每一份3個點做3點(2階)拉格朗日插值
% 優(yōu)勢: 避免高階牛頓-科茨(高階)出現(xiàn)龍貝格現(xiàn)象, 因此等分越多越好!
% 以f(x) = 4/(1+x^2)函數(shù)為例

clear; clc;

syms x;
f = 4/(1+x^2);   % 每次這里修改不同函數(shù)即可

up = double(input('輸入積分上限:'));
low = double(input('輸入積分下限:'));
n = double(input('將積分區(qū)間n等分(偶數(shù)):'));
m = n/2;
% 真實定積分的結(jié)果: 
R = int(f,x,low,up); 

fmiddle1 = zeros(1,m);    % 對大括號中第一項的記錄
fmiddle2 = zeros(1,m-1);  % 對大括號中第二項的記錄
for num1 = 1:m
    x = low + (2*num1-1)*(up-low)/(2*m);
    fmiddle1(num1) = double(subs(f));
end

for num2 = 1:m-1
    x = low + (2*num2)*(up-low)/(2*m);
    fmiddle2(num2) = double(subs(f));
end

x = up;
fup = double(subs(f));
x = low;
flow = double(subs(f));
% 復化梯形公式近似結(jié)果:
FS = (up-low)/(6*m)*( flow + 4*sum(fmiddle1) + 2*sum(fmiddle2) + fup ); 

fprintf('復化梯形求積公式近似結(jié)果為: %.7f\n',FS);
fprintf('真實結(jié)果為: %f\n',R);

（3）復化梯形公式加密：名字看上去很復雜扒最，其實就是一個快速計算復化梯形公式當前分區(qū)數(shù)加倍后的新復化梯形公式結(jié)果的公式。就是一個計算規(guī)律而已华嘹。具體含義如下：

當區(qū)間n等分時吧趣，復化梯形結(jié)果為：

$T(h)$

現(xiàn)將每個小區(qū)間對分，即將原區(qū)間劃分成為2n份耙厚，新分區(qū)下每個區(qū)間長度為：h = (b-a)/(2n)强挫；此時新的結(jié)果有一個快速計算的方法(注意：式中紅色n是原始的等分數(shù)，常數(shù)不會變薛躬；h每次都變)：

$T(\frac{h}{2}) = \frac{T(h)}{2} + \frac{h}{2}\sum_{i=1}^{\color{red}{n}}f(a + (2i-1)\frac{h}{2})$

現(xiàn)將新的每個小區(qū)間再對分俯渤，即將原區(qū)間劃分成為4n份，新分區(qū)下每個區(qū)間長度為：h = (b-a)/(4n)型宝；此時新的結(jié)果為：

$T(\frac{h}{2^2}) = \frac{T(\frac{h}{2})}{2} + \frac{h}{2^2}\sum_{i=1}^{\color{red}{2n}}f(a + (2i-1)\frac{h}{2^2})$

$\vdots$

通式為：

$T(\frac{h}{2^{k+1}}) = \frac{T(\frac{h}{2^k})}{2} + \frac{h}{2^{k+1}}\sum_{i=1}^{\color{red}{2^k}}f(a + (2i-1)\frac{h}{2^{k+1}}) \quad (k=0,1,2,\cdots)$

總結(jié)一句：復化梯形加密八匠，其實就是基于現(xiàn)有分區(qū)情況下的復化梯形結(jié)果，快速推算每個小區(qū)間再對分的新結(jié)果趴酣；新結(jié)果相當于再原基礎上又增加了2倍的插值節(jié)點梨树。就是一個快推的規(guī)律。

% 復化梯形公式加密
% 作用: 快速計算"等分"后再在"每個小份不斷對分"的加速算法! 就是為了加速計算的
% 方式: 知道1等分后, 用加密方法可以快速計算任意2^k等分的結(jié)果
% 以f(x) = 4/(1+x^2)函數(shù)為例

clear; clc

syms x;
f = 4/(1+x^2);   % 每次這里修改不同函數(shù)即可

up = double(input('輸入積分上限:'));
low = double(input('輸入積分下限:'));
k = double(input('輸入目標對分/加密次數(shù):'));
% 真實定積分的結(jié)果: 
R = int(f,x,low,up);  

% 先計算1等分的結(jié)果, 后面加密高等分都是基于此推出的岖寞！
x = up;
fup = double(subs(f));
x = low;
flow = double(subs(f));
Torg = (up - low)/2*(fup+flow);  % 0等分(梯形公式結(jié)果)

% 復化梯形加密:
T = zeros(1,k+1);  % 記錄每一次加密的結(jié)果
T(1) = Torg;       % 把1等分的放到第一個, 方便后面的循環(huán)
h = up - low;  % 區(qū)間原長
tmp = 1;       % 記錄T索引用的
for numk = 0:k
    fmiddle = zeros(1,2^(numk));  
    % 小循環(huán)里所有和k有關(guān)都是變化的!都和當前numk有關(guān):
    for numi = 1:2^(numk)
        x = low + (2*numi-1)*h/2^(numk+1);
        fmiddle(numi) = double(subs(f));
    end
    T(tmp+1) = 0.5*T(tmp) + h/2^(numk+1)*sum(fmiddle);
    tmp = tmp + 1;
    if tmp == k+1    % 已經(jīng)滿了,可以結(jié)束了
        break;
    end
end

% 加密過程中的一系列結(jié)果都在T中, 打印只顯示最后精度最高的那個T(k+1):
fprintf('復化梯形加密為%d等分后,當前近似結(jié)果為: %.7f\n',2^k,T(k+1));
fprintf('真實結(jié)果為: %.10f\n',R);

（4）龍貝格公式：就是在復化梯形加密之后抡四，再做一個加速操作并進一步提高精度！加速的內(nèi)涵其實是改變多項式相加時的各個系數(shù)慎璧。一般設 $T_0(\frac{h}{2^k})$ 表示k次加密后復化梯形公式的值(T的下標是0即表示沒有加速)床嫌；以 $T_m(\frac{h}{2^k})$ 表示在k加密的基礎上又做了m次加速后的結(jié)果跨释。加密與加速的遞推公式為：

$T_m(\frac{h}{2^k}) = \frac{4^mT_{\color{red}{m-1}}(\frac{h}{2^{k+1}})-T_\color{red}{m-1}(\frac{h}{2^k})}{4^m-1}$

其中加密次數(shù)k和加速次數(shù)m的數(shù)值范圍為：

$k = 0,1,2,\cdots \quad \quad m = 1,2,\cdots,k+1$

過程舉例：先依次做0,1,2,3,4次加密(都做胸私！)，然后在這4個加密結(jié)果基礎上鳖谈，依次做1,2,3,4次加速岁疼；最后得到的在4次加密基礎上的4次加速，精度是最高的。先做完所有的加密操作捷绒，再做所有的加速操作瑰排，兩個過程是先后獨立進行的。

% 龍貝格公式——在加密復化梯形公式的基礎上再循環(huán)加速(提高精度)!
% 作用: 在加密復化梯形公式基礎上, 再提高精度
% 方式: 先全部加密完后, 再依次進行加速——兩個過程是"獨立"進行的暖侨！
% 以f(x) = 4/(1+x^2)函數(shù)為例

clear; clc
format long;

syms x;
f = 4/(1+x^2);   % 每次這里修改不同函數(shù)即可

up = double(input('輸入積分上限:'));
low = double(input('輸入積分下限:'));
k = double(input('輸入目標對分/加密次數(shù):'));
% 真實定積分的結(jié)果: 
R = int(f,x,low,up);  

% 先計算1等分的結(jié)果, 后面加密高等分都是基于此推出的椭住！
x = up;
fup = double(subs(f));
x = low;
flow = double(subs(f));
Torg = (up - low)/2*(fup+flow);  % 0等分(梯形公式結(jié)果)

% 復化梯形加密:
T = zeros(1,k+1);  % 記錄每一次加密的結(jié)果
T(1) = Torg;       % 把1等分的放到第一個, 方便后面的循環(huán)
h = up - low;  % 區(qū)間原長
tmp = 1;       % 記錄T索引用的
for numk = 0:k
    fmiddle = zeros(1,2^(numk));  
    % 小循環(huán)里所有和k有關(guān)都是變化的!都和當前numk有關(guān):
    for numi = 1:2^(numk)
        x = low + (2*numi-1)*h/2^(numk+1);
        fmiddle(numi) = double(subs(f));
    end
    T(tmp+1) = 0.5*T(tmp) + h/2^(numk+1)*sum(fmiddle);
    tmp = tmp + 1;
    if tmp == k+1    % 已經(jīng)滿了,可以結(jié)束了
        break;
    end
end

% 龍貝格/理查森加速: 用表格/矩陣表示更好編程
Rbg = zeros(k+1,k+1);
Rbg(:,1) = T;  % 第一列初始化(為加速前,即所有的加密點結(jié)果)

m = 1;
for col = 2:k+1
    for row = 1:k+1-m
        Rbg(row,col) = ( 4^m * Rbg(row+1,col-1) - Rbg(row,col-1) )/( 4^m-1 );
    end
    m = m + 1;
end

% 加密過程中的一系列結(jié)果都在T中, 打印只顯示最后精度最高的那個T(k+1):
fprintf('復化梯形加密為%d等分后,當前近似結(jié)果為: %.10f\n',2^k,T(k+1));
fprintf('基于當前加密后的龍貝格加速結(jié)果為:%.10f\n',Rbg(1,k+1));
fprintf('真實結(jié)果為: %.10f\n',R);

三、高斯型積分公式

前文我們提到：拉格朗日型求積公式要想提高精度字逗，只能不斷增加插值節(jié)點京郑！它存在的一個問題是：選擇的插值點"效率不高"，即都是一些沒有實際內(nèi)涵和意義的點(比如均分點葫掉、等比分點等)些举。可以想象到：如果插值節(jié)點都是"精心挑選"的俭厚，很能代表原函數(shù)特征和變化趨勢的有意義的點户魏，那對這些點再做n階拉格朗日插值精度肯定會再次提高！高斯型積分公式挪挤，就是為了精心挑選插值節(jié)點叼丑。

拉格朗日型：大量增加效率不高插值節(jié)點；然后非重疊2點間做線性拉格朗日插值扛门；
高斯型：插值節(jié)點數(shù)目一定時幢码，通過精心調(diào)整節(jié)點的位置，提高每個節(jié)點的效率尖飞；然后做n+1個節(jié)點的n階拉格朗日插值症副。

高斯型缺點：用的還是高階拉格朗日插值！Ｕ贞铣！本文還是采用這種方式。
高斯型改進：在選擇了好的插值點后沮明，再使用分段線性插值辕坝。

高斯型求積公式要用到正交多項式相關(guān)知識。簡言之就是求n階正交多項式的零點(插值節(jié)點 $x_i$ )和系數(shù)(求積系數(shù) $A_i$ )荐健；節(jié)點和系數(shù)的計算公式如下：

$P_n(x_i) = 0$

$A_i = \int_{a}^酱畅p(x)l_i(x_i)dx \quad\quad (i=0,1,2,\cdots,n)$

節(jié)點求解：另n階正交多項式為0，求這個方程的零點(這就是"精心挑選"的過程)江场。幾階正交多項式就會用幾個零點纺酸。
系數(shù)求解： $p(x)$ 是選擇的正交多項式序列所對應的權(quán)函數(shù)； $l_i(x)$ 是n階拉格朗日插值多項式(插值節(jié)點就是上面求得的 $x_i$ )址否。

有了插值節(jié)點和對應的系數(shù)后餐蔬，原定積分式可以進行如下近似計算：

$\int_{a}^p(x)f(x)dx ≈ \sum_{i=0}^{n}A_if(x_i)$

注意：左邊待積函數(shù)中必須帶有所選的正交多項式對應的權(quán)函數(shù)！樊诺！個人覺得這樣雖然限制了高斯型積分的普遍性仗考，但是特殊問題可以選用特殊的正交多項式序列和對應權(quán)函數(shù)來解決，其實影響不大词爬。注意到：如果權(quán)函數(shù) $p(x)=1$ 的話秃嗜，那么相當于任何原待積函數(shù)都有這個權(quán)函數(shù)，也就是說該權(quán)函數(shù)所對應的那種正交多項式序列最具有通用性顿膨！這就是我最推薦使用的：高斯-勒讓德求積公式痪寻。高斯-勒讓德多項式在這里查看。

高斯-勒讓德節(jié)點虽惭、系數(shù)橡类、求積一套流程的Matlab程序如下：

% 任意個數(shù)插值節(jié)點的高斯-勒讓德求積公式:
% 思路: 其內(nèi)涵還是n階拉格朗日插值! 只不過原先的拉格朗日插值點都是等間距均勻分布的那種,
%       高斯-勒讓德就是點數(shù)一定時, 挑選那些能讓精度提高的點(選最好的點)!而不再使用無特殊意義的等間距分隔點。
%       補充: 高斯-勒讓德"不能"與前面的方法混用(比如選最好點后再龍貝格)!
%             因為前面所有的方法要求"必須是等間距"的點才能寫出"那樣的公式"!
% 優(yōu)點: 小區(qū)間/積分限上,限制了只能有幾個控制點時, 用高斯積分效果不錯! 
% 缺點: 大區(qū)間上, 高斯-勒讓德的求積節(jié)點必須得增多! 因為內(nèi)涵還是n階拉格朗日插值!——這個缺點非常大!
% 建議: 沒有控制點限制時, 我還是推薦用龍貝格! 因為它是萬能適用的高精度! 

clear; clc;
format long;

syms x;
n = double(input('輸入使用幾點(n)的高斯-勒讓德插值:'));
% n點插值的高斯-勒讓德多項式P:
f = x^2 -1;
fprintf('%d點高斯-勒讓德多項式為:\n',n)
P = vpa(1/(2^n*factorial(n)) * diff(f^n,x,n))
% n點高斯-勒讓德插值節(jié)點Xi:
Xi = sort(double(solve(P)))';

% n點高斯-勒讓德插值節(jié)點對應的插值系數(shù)Ai:
xnum = Xi;
% l用來記錄拉格朗日多項式的: n點拉格朗日插值多項式有n個基函數(shù)! 
% 高斯-勒讓德求積系數(shù)就是對"每個多項式/基函數(shù)"求其[-1 1]定積分——
% 插值點數(shù)芽唇、多項式/基函數(shù)個數(shù)顾画、求積系數(shù)個數(shù)一樣，相互對應匆笤。
l = sym(zeros(1,n));  
Ai = zeros(1,n);
for m = 1:n
    l(m) = prod(x - xnum([1:m-1 m+1:n]))/prod(xnum(m) - xnum([1:m-1 m+1:n])); 
    Ai(m) = double(int(l(m),x,-1,1));
end

fprintf('%d點高斯-勒讓德插值節(jié)點為:\n',n);
Xi
fprintf('%d點高斯-勒讓德插值節(jié)點對應的系數(shù)為:\n',n);
Ai


% 對具體函數(shù)的積分近似: 
syms t;
up = double(input('輸入積分上限:'));
low = double(input('輸入積分下限:'));
% 積分限x必須是[-1 1], 因此這里統(tǒng)一做一下轉(zhuǎn)換, 即任何限都可以:
t = (up-low)/2*x + (up+low)/2;
y = 4/(1+t^2);  % 每次修改這里的函數(shù)即可, 注意這里自變量是t (修改√√√)

Result = 0;  % 記錄最終近似結(jié)果
for num = 1:n
    x = Xi(num);
    y_tmp = (up-low)/2*double(subs(y));
    Result = Result + Ai(num)*y_tmp;
end

syms w;
y_r = 4/(1+w^2);  % 實際的結(jié)果, 注意這里自變量是w  (修改√√√)
R = int(y_r,w,low,up);
fprintf('%d點高斯-勒讓德積分近似結(jié)果為:%.10f\n', n, Result);
fprintf('真實結(jié)果為:%.10f\n', R);

需要注意一點的是：高斯-勒讓德的積分限必須是 $[-1,1]$ 研侣！如果你給的積分限不是這個，程序中會自動幫你轉(zhuǎn)化到這個范圍當中炮捧。

3.1 其他高斯型求積公式

其他常用的還有：1. 高斯-拉蓋爾求積庶诡；2. 高斯-埃爾米特求積。多項式參考這里咆课。

（1）高斯-拉蓋爾求積：積分區(qū)間為 $[0,+\infty)$ 末誓，權(quán)函數(shù)為：

$p(x) = e^{-x}$

使用時：原被函數(shù)函數(shù)中必須有一項相乘的 $e^{-x}$ 才能使用這種。

% 任意個數(shù)插值節(jié)點的高斯-拉蓋爾求積公式:
% 說明等信息同高斯-勒讓德
% 用處: 原始待積函數(shù)形式中必須有一個exp(-x)項! 即f(x) = exp(-x)*g(x)
% 注意: 近似計算時用的是g(x)! 即Ai*g(xi)! 前面那個exp(-x)在近似計算中是不用的!

clear; clc;
format long;

syms x;
n = double(input('輸入使用幾點(n)的高斯-拉蓋爾插值:'));
% n點插值的高斯-拉蓋爾多項式P:
f = exp(-x)*x^n;
fprintf('%d點高斯-拉蓋爾多項式為:\n',n)
L = vpa( exp(x) * diff(f,x,n) )
% n點高斯-拉蓋爾插值節(jié)點Xi:
Xi = sort(double(solve(L)))';

% n點高斯-拉蓋爾插值節(jié)點對應的插值系數(shù)Ai:
xnum = Xi;
% l用來記錄拉格朗日多項式的: n點拉格朗日插值多項式有n個基函數(shù)! 
% 高斯-拉蓋爾求積系數(shù)就是對"每個多項式/基函數(shù)"求其[0 +inf]定積分——
% 插值點數(shù)书蚪、多項式/基函數(shù)個數(shù)喇澡、求積系數(shù)個數(shù)一樣，相互對應殊校。
l = sym(zeros(1,n));  
Ai = zeros(1,n);
for m = 1:n
    l(m) = prod(x - xnum([1:m-1 m+1:n]))/prod(xnum(m) - xnum([1:m-1 m+1:n])); 
    l_tmp = exp(-x)*l(m);
    Ai(m) = double(int(l_tmp,x,0,+inf));
end

fprintf('%d點高斯-拉蓋爾插值節(jié)點為:\n',n);
Xi
fprintf('%d點高斯-拉蓋爾插值節(jié)點對應的系數(shù)為:\n',n);
Ai


% 對具體函數(shù)的積分近似: 高斯-拉蓋爾積分限很簡單,必須是[0,+inf]
syms t;
gt = sin(t);    % 每次修改這里的函數(shù)即可, 注意exp(-t)這一項必須有
y = exp(-t)*gt;   
up = +inf;
low = 0;
R = int(y,t,low,up);  % 實際結(jié)果

Result = 0;  % 記錄最終近似結(jié)果
for num = 1:n
    t = Xi(num);
    y_tmp = double(subs(gt));
    Result = Result + Ai(num)*y_tmp;
end

fprintf('%d點高斯-拉蓋爾積分近似結(jié)果為:%.10f\n', n, Result);
fprintf('真實結(jié)果為:%.10f\n', R);

（2）高斯-埃爾米特求積：積分區(qū)間為 $(-\infty,+\infty)$ 晴玖，權(quán)函數(shù)為：

$p(x) = e^{x^2}$

使用時：原被函數(shù)函數(shù)中必須有一項相乘的 $e^{x^2}$ 才能使用這種。

% 任意個數(shù)插值節(jié)點的高斯-埃爾米特求積公式:
% 說明等信息同高斯-勒讓德
% 用處: 原始待積函數(shù)形式中必須有一個exp(-x^2)項! 即f(x) = exp(-x^2)*g(x)
% 注意1: 近似計算時用的是g(x)! 即Ai*g(xi)! 前面那個exp(-x^2)在近似計算中是不用的!
% 注意2: 埃爾米特特殊! 當插值節(jié)點大于4個后, 只能用"奇數(shù)個數(shù)"的插值點!!!!!!!!
%        因為"解插值節(jié)點的方程"的解此時個數(shù) < n 也就不能再往下計算了为流。

clear; clc;
format long;

syms x;
n = double(input('輸入使用幾點(n)的高斯-埃爾米特插值:'));
% n點插值的高斯-埃爾米特多項式P:
f = exp(-x^2);
fprintf('%d點高斯-埃爾米特多項式為:\n',n)
H = vpa( (-1)^(n-1) * exp(x^2) * diff(f,x,n) )
% n點高斯-埃爾米特插值節(jié)點Xi:
Xi = sort(double(solve(H)))'

% n點高斯-埃爾米特插值節(jié)點對應的插值系數(shù)Ai:
xnum = Xi;
% l用來記錄拉格朗日多項式的: n點拉格朗日插值多項式有n個基函數(shù)! 
% 高斯-埃爾米特求積系數(shù)就是對"每個多項式/基函數(shù)"求其[-inf +inf]定積分——
% 插值點數(shù)呕屎、多項式/基函數(shù)個數(shù)、求積系數(shù)個數(shù)一樣敬察，相互對應秀睛。
l = sym(zeros(1,n));  
Ai = zeros(1,n);
for m = 1:n
    l(m) = prod(x - xnum([1:m-1 m+1:n]))/prod(xnum(m) - xnum([1:m-1 m+1:n])); 
    l_tmp = exp(-x^2)*l(m);
    Ai(m) = double(int(l_tmp,x,-inf,+inf));
end

fprintf('%d點高斯-埃爾米特插值節(jié)點為:\n',n);
Xi
fprintf('%d點高斯-埃爾米特插值節(jié)點對應的系數(shù)為:\n',n);
Ai


% 對具體函數(shù)的積分近似: 高斯-拉蓋爾積分限很簡單,必須是[-inf,+inf]
syms t;
gt = sin(t)^2;    % 每次修改這里的函數(shù)即可, 注意exp(-t)這一項必須有
y = exp(-t^2)*gt;   
up = +inf;
low = -inf;
R = int(y,t,low,up);  % 實際結(jié)果

Result = 0;  % 記錄最終近似結(jié)果
for num = 1:n
    t = Xi(num);
    y_tmp = double(subs(gt));
    Result = Result + Ai(num)*y_tmp;
end

fprintf('%d點高斯-埃爾米特積分近似結(jié)果為:%.10f\n', n, Result);
fprintf('真實結(jié)果為:%.10f\n', R);

注意：高斯-埃爾米特的使用，當插值節(jié)點大于4個后静汤，只能用奇數(shù)個插值節(jié)點琅催。

最后編輯于：2019.04.07 10:36:09

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個濱河市虫给，隨后出現(xiàn)的幾起案子藤抡，更是在濱河造成了極大的恐慌，老刑警劉巖抹估，帶你破解...
沈念sama閱讀 206,482評論 6贊 481
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件缠黍，死亡現(xiàn)場離奇詭異，居然都是意外死亡药蜻，警方通過查閱死者的電腦和手機瓷式，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 88,377評論 2贊 382
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來语泽，“玉大人贸典，你說我怎么就攤上這事□饴眩” “怎么了廊驼？”我有些...
開封第一講書人閱讀 152,762評論 0贊 342
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長惋砂。經(jīng)常有香客問我妒挎，道長，這世上最難降的妖魔是什么西饵？我笑而不...
開封第一講書人閱讀 55,273評論 1贊 279
?港島之戀（遺憾婚禮）
正文為了忘掉前任酝掩，我火速辦了婚禮，結(jié)果婚禮上眷柔，老公的妹妹穿的比我還像新娘期虾。我一直安慰自己，他們只是感情好驯嘱，可當我...
茶點故事閱讀 64,289評論 5贊 373
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布彻消。她就那樣靜靜地躺著，像睡著了一般宙拉。火紅的嫁衣襯著肌膚如雪宾尚。梳的紋絲不亂的頭發(fā)上，一...
開封第一講書人閱讀 49,046評論 1贊 285
城市分裂傳說
那天谢澈，我揣著相機與錄音煌贴，去河邊找鬼。笑死锥忿，一個胖子當著我的面吹牛牛郑，可吹牛的內(nèi)容都是我干的。我是一名探鬼主播敬鬓，決...
沈念sama閱讀 38,351評論 3贊 400
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼淹朋，長吁一口氣：“原來是場噩夢啊……” “哼笙各！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起础芍，我...
開封第一講書人閱讀 36,988評論 0贊 259
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤杈抢，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后仑性，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體惶楼，經(jīng)...
沈念sama閱讀 43,476評論 1贊 300
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 35,948評論 2贊 324
?白月光啟示錄
正文我和宋清朗相戀三年诊杆，在試婚紗的時候發(fā)現(xiàn)自己被綠了歼捐。大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
茶點故事閱讀 38,064評論 1贊 333
活死人
序言：一個原本活蹦亂跳的男人離奇死亡晨汹，死狀恐怖豹储，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情淘这，我是刑警寧澤颂翼，帶...
沈念sama閱讀 33,712評論 4贊 323
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站慨灭，受9級特大地震影響朦乏，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜氧骤，卻給世界環(huán)境...
茶點故事閱讀 39,261評論 3贊 307
男人毒藥：我在死后第九天來索命
文/蒙蒙一呻疹、第九天我趴在偏房一處隱蔽的房頂上張望。院中可真熱鬧筹陵，春花似錦刽锤、人聲如沸。這莊子的主人今日做“春日...
開封第一講書人閱讀 30,264評論 0贊 19
一樁弒父案并思，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽。三九已至语稠，卻和暖如春宋彼，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背仙畦。一陣腳步聲響...
開封第一講書人閱讀 31,486評論 1贊 262
情欲美人皮
我被黑心中介騙來泰國打工输涕，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留，地道東北人慨畸。一個月前我還...
沈念sama閱讀 45,511評論 2贊 354
代替公主和親
正文我出身青樓莱坎，卻偏偏與公主長得像，于是被迫代替她去往敵國和親寸士。傳聞我的和親對象是個殘疾皇子檐什，可洞房花燭夜當晚...
茶點故事閱讀 42,802評論 2贊 345