对于积分,只要找到被积函数
的原函数
,便有下列牛顿--莱布尼兹公式:
但实际上采用这种方法往往有困难,一是有些原函数不能用初等函数表达,有些计算过于麻烦。因为我们可以用数值的方法计算。
积分中值定理告诉我们,在区间 [a,b]内存在一点 ,使得
.
就是说,底为b-a而高为的矩形面积恰等于所求曲边梯形的面积 I ,问题在于点
一般是不知道的。但我们可以将
称为区间上的平均高度,这样我们只要对平均高度提供一种算法,便可以获得一种数值求积方法。
一.小区间上的数值积分
1.梯形公式
我们找一个函数,即最简单的一次函数 来拟合
,用两端点 f(a),f(b) 的算术平均值作为平均高度的近似值,那么
,那么求积公式为
,这就是梯形公式。
代码如下:
%f(x)=2*x
a=1;b=1.1;
disp("梯形公式解");
t=((b-a)/2)*(hanshu(a)+hanshu(b))
disp("MATLAB自带函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)
function[f]=hanshu(x)
f=x*x;
end
结果:
梯形公式的代数精度为1,不再赘述。
2.辛普森公式
如果我们不用最简单的函数,而改用二次函数,那么便是辛普森公式,直接写公式
disp("sinpson解")
t=((b-a)/6)*(hanshu(a)+hanshu(b)+4*hanshu((a+b)/2))
disp("函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x)
f=x*x;
end
结果:
同样的函数,辛普森公式的解就比梯形公式的解更精确了,当然计算复杂度也高了。
二:一般的数值积分
1. 复合的梯形公式
一般的区间主要是切分成小区间,再利用梯形公式,
对区间 [a,b]做切分,分成m块,Xi=(b-a)/m;
分成m块:
利用梯形公式:(中间的边需要算2遍,所以乘2)
代码如下:
disp("复合的梯形公式,分成2块")
m=2;step=(b-a)/m;
t3=((b-a)/(2*m))*(hanshu(a)+hanshu(b)+xun(a,step,b))
disp("复合的梯形公式,分成5块")
m=5;step=(b-a)/m;
t3=((b-a)/(2*m))*(hanshu(a)+hanshu(b)+xun(a,step,b))
disp("函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x)
f=x*x;
end
function[t]=xun(a,step,b)
t=0;
for i=a+step:step:b-step
t=t+hanshu(i);
end
t=2*t;
end
结果:
2. 复合的辛普森公式
直接写公式:
代码如下:
disp("复合的sinpson公式,分成1块")
m=1;step=(b-a)/m;
t3=((b-a)/(6*m))*(hanshu(a)+hanshu(b)+xun(a,step,b)+xu(a,step,b))
disp("复合的sinpson公式,分成2块")
m=2;step=(b-a)/m;
t3=((b-a)/(6*m))*(hanshu(a)+hanshu(b)+xun(a,step,b)+xu(a,step,b))
disp("函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x)
f=x*x;
end
function[t]=xun(a,step,b)
t=0;
for i=a+step:step:b-step
t=t+hanshu(i);
end
t=2*t;
end
function[xx]=xu(a,step,b)
xx=0;
for i=a+step/2:step:b
xx=xx+hanshu(i);
end
xx=xx*4;
end
结果如下:
----------------------------分割线----------------------------------------
加了注释的代码:
%梯形公式的代码
clc;clear all;
disp("被积函数: y=x*x");
disp("梯形公式解");
a=input('积分下限:a=');
b=input('积分上限:b=');
ft1=((b-a)/2)*(hanshu(a)+hanshu(b))
disp("MATLAB自带函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
%复合的梯形公式
disp("复合的梯形公式,分成多少块")
m=input('m=');
step=(b-a)/m; % m为分成几份,step为梯形的高
ft2=((b-a)/(2*m))*(hanshu(a)+hanshu(b)+xun(a,step,b))
disp("MATLAB自带函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x)%被积函数y
f=x*x;
end
function[t]=xun(a,step,b) %求复合梯形公式内部梯形被求2次的上底和下底
t=0;
for i=a+step:step:b-step %循环每次加一个步长
t=t+hanshu(i);
end
t=2*t;
end
%sinpson公式的代码
clc;clear all;
disp("sinpson公式解");
a=input('积分下限:a=');
b=input('积分上限:b=');
ft1=((b-a)/6)*(hanshu(a)+hanshu(b)+4*hanshu((a+b)/2))
disp("MATLAB自带函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
disp("复合的sinpson公式,分成多少块")
m=input('m=');
step=(b-a)/m;
ft2=((b-a)/(6*m))*(hanshu(a)+hanshu(b)+xun(a,step,b)+xu(a,step,b))
disp("MATLAB自带函数解")
syms x;
y=x*x;
x=int(y,x,a,b);
x=vpa(x)%转小数
function[f]=hanshu(x) %被积函数
f=x*x;
end
function[t]=xun(a,step,b) %梯形重复计算的上底和下底
t=0;
for i=a+step:step:b-step
t=t+hanshu(i);
end
t=2*t;
end
function[xx]=xu(a,step,b) %sinpson公式中的 4*f((a+b)/2)
xx=0;
for i=a+step/2:step:b
xx=xx+hanshu(i);
end
xx=xx*4;
end