8 数值积分与数值微分

8 数值积分与数值微分

8.1 例题解答

例 8.1 给定积分解:

先输入主要初始参数

>>a=0.5; >>b=1;

>>f=inline('x^(1/2)');

%梯形公式

>>I1=(b-a)/2*(feval(f,a)+feval(f,b)) I1 =

0.426776695296637

%simpson公式

>>I2=(b-a)/6*(feval(f,a)+4*feval(f,(a+b)/2)+feval(f,b)) I2 =

0.430934033027025 %Cotes公式(n=4) >>tc=0;

>>C0=[7 32 12 32 7]; >>for i=0:4

tc=tc+C0(i+1)*feval(f,a+i*(b-a)/4); end

>>I3=(b-a)/90*tc I3 =

0.430964070495876

%准确值

>>I=int(char(f),a,b) >>vpa(I) I =

-1/6*2^(1/2)+2/3 ans =

0.43096440627115082519971854596505

?10.5xdx,分别用梯形公式、Simpson公式、Cote公式作近似计算.

例 8.2 对积分I?sinx?4?0xdx,为使其精度达到10.

1若用复化梯形公式,应将[0,1]多少等分?

若用复化Simpson公式,应将[0,1]多少等分? 解:

直接按余项计算即可. 复化梯形公式的余项为:

En(f)??复化Simpson公式余项为:

b?a2''hf(?),12??(a,b)

b?a?h?(4)En(f)????f(?),180?2?对于f(x)?4??(a,b),f(x)?C4[a,b]

sinx,在课本中我们已证得以下不等式成立: x1f(k)(x)?

k?1直接利用上述不等式关系解答本题. 先输入误差精度: >>Eps=1E-4

Eps =

1.000000000000000e-004

(1) 复化梯形公式

>>h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1))) %先求出步长

h1 =

0.060000000000000

>>N1=ceil(1/h1) %向上取整,得到等分区间数 N1 = 17

故可将区间17等分即可达到所要求的精度. (2) 复化Simpson公式

>>h2=power(Eps/abs(-(1-0)/180*1/(1+4)),1/4) %先求出步长 h2 =

0.547722557505166

>>N2=ceil(1/h2) %向上取整,得到等分区间数 N2 = 2

故可将区间2等分即可达到所要求的精度.

? 扩展:

1) Matlab中复化梯形公式命令为I=trapz(x,y),复化Simpson公式命令为quad().

2) Matlab中有四个取整函数,分别为ceil(),floor(),fix(),round(),分别表示向正无穷大方向取

整、向负无穷大方向取整、向靠近零方向舍入和四舍五入.

例 8.3 对积分I?sinx?6dx,利用变步长方法求其近似值,使其精度达到. ??10?0x1解：

利用变步长法前先建立三种变步长复化积分公式的函数. 注意在Matlab中直接用sin(0)/0得不到1，此函数名为limit().

先建立三种复化公式的函数文件,它们分别为复化梯形公式trap.m、复化Simpson公式为simpson.m、Cotes公式为cotes.m,三个函数的源程序如下: (1) 复化梯形公式trap.m

function T=trap(f,a,b,n) %trap.m

%复化梯形公式求积分值 %f为积分函数 %[a,b]为积分区间 %n是等分区间份数

h=(b-a)/n;%步长 T=0;

for k=1:(n-1) x0=a+h*k; T=T+limit(f,x0); end

T=h*(limit(f,a)+limit(f,b))/2+h*T; T=double(T);

(2) 复化Simpson公式simpson.m:

function S=simpson(f,a,b,n) %simpson.m

%Simpson公式求积分值 %f为积分函数 %[a,b]为积分区间 %n是等分区间份数

h=(b-a)/(2*n);%步长 s1=0; s2=0; for k=1:n

sin(0)0matlab?NAN，因此解此题时我们改用求极限的方法得到函数值,

x0=a+h*(2*k-1); s1=s1+limit(f,x0); end

for k=1:(n-1)

x0=a+h*2*k; s2=s2+limit(f,x0); end

S=h*(limit(f,a)+limit(f,b)+4*s1+2*s2)/3; S=double(S);

(3) 复化Cotes公式cotes.m:

function C=cote(f,a,b,n) %cote.m

%复化cotes公式求积分值 %f为积分函数 %[a,b]为积分区间 %n是等分区间份数

h=(b-a)/n;%步长 C=0;

for i=1:(n-1) x0=a+i*h;

C=C+14*limit(f,x0); end

for k=0:(n-1) x0=a+h*k;

s=32*limit(f,x0+h*1/4)+12*limit(f,x0+h*1/2)+32*limit(f,x0+h*3/4); C=C+s; end

C=C+7*(limit(f,a)+limit(f,b)); C=C*h/90; C=double(C);

再编写主程序调用这三个函数,主程序名为ex8_3.m,源程序如下:

%ex8_3.m clc; syms x;

f=sym('sin(x)/x');

a=0;b=1;%积分上下限

n=20;%作1,2,3,…,20次区间等分 %复化梯形公式 T=zeros(n,1);

for i=1:n

T(i)=trap(f,a,b,i); end

%复化Simpson公式; S=zeros(n,1); for i=1:n

S(i)=simpson(f,a,b,i); end

%复化Cotes公式 C=zeros(n,1); for i=1:n

C(i)=cote(f,a,b,i); end %准确值 I=int(f,a,b); I=double(I);

%画图作出直观观察 x=[]; x=1:n; figure;

plot(x,ones(1,n)*I,'-'); hold on;

plot(x,T','r--','LineWidth',2); plot(x,S','m.-','LineWidth',1); plot(x,C','c:','LineWidth',1.5); grid on;

title('三种复化公式积分效果对比图');

legend('准确值曲线','复化梯形公式','复化Simpson公式','复化Cotes公式'); hold off;

disp(' 复化梯形公式 复化Simpson公式 复化Cotes公式'); disp([T S C]);

在Matlab命令窗口中输入ex8_3,得到如下积分结果:

复化梯形公式 复化Simpson公式 复化Cotes公式 0.920735492403948 0.946145882273587 0.946083004063674 0.939793284806177 0.946086933951794 0.946083069350917 0.943291429132337 0.946083831311699 0.946083070278278 0.944513521665390 0.946083310888472 0.946083070351379 0.945078780953402 0.946083168838073 0.946083070363043 0.945385730766859 0.946083117842867 0.946083070365797 0.945570776256246 0.946083095989403 0.946083070366633 0.945690863582701 0.946083085384948 0.946083070366936

0.945773188549752 0.946083079742053 0.946083070367061 0.945832071866905 0.946083076517732 0.946083070367118 0.945875637119198 0.946083074567937 0.946083070367147 0.945908771073865 0.946083073333114 0.946083070367161 0.945934556488475 0.946083072520476 0.946083070367170 0.945955016056114 0.946083071968057 0.946083070367174 0.945971521552985 0.946083071581964 0.946083070367177 0.945985029934386 0.946083071305562 0.946083070367179 0.945996225242376 0.946083071103489 0.946083070367180 0.946005606943978 0.946083070952998 0.946083070367181 0.946013546623966 0.946083070839065 0.946083070367182 0.946020325355025 0.946083070751532 0.946083070367182

三种复化公式积分效果对比图0.955准确值曲线复化梯形公式复化Simpson公式复化Cotes公式 0.950.9450.940.9350.930.9250.92 02468101214161820 由以上结果可看到复化梯形公式有一个上升接近准确值过程,而复化Simpson公式积分结果和复化Cotes公式积分的结果基本上和准确值的曲线重叠在一块,可见它们的精度是相当高的.

例 8.4 用Romberg积分法计算I??解：

编写Romberg积分法的函数M文件,源程序如下(romberg.m):

function [I,T]=romberg(f,a,b,n,Eps) %Romberg积分计算 %f为积分函数 %[a,b]为积分区间

sinxdx,精度??10?6. 0x1%n+1是T数表的列数目 %Eps为迭代精度

%返回值中I为积分结果,T是积分表

if nargin<5 Eps=1E-6; end m=1; h=(b-a); err=1; j=0;

T=zeros(4,4);

T(1,1)=h*(limit(f,a)+limit(f,b))/2; while ((err>Eps) & (j

for p=1:m

x0=a+h*(2*p-1); s=s+limit(f,x0); end

T(j+1,1)=T(j,1)/2+h*s; m=2*m; for k=1:j

T(j+1,k+1)=T(j+1,k)+(T(j+1,k)-T(j,k))/(4^k-1); end

err=abs(T(j,j)-T(j+1,k+1)); end

I=T(j+1,j+1); if nargout==1 T=[]; end

将上述源程序另存为romberg.m后,进入计算:

>>syms x;%创建符号变量

>>f=sym('sin(x)/x') %符号函数 f = sin(x)/x

>> [I,T]=romberg(f,0,1,3,1E-6) %积分计算I = 0.9461 T =

0.9207 0 0 0 0

0.9398 0.9461 0 0 0 0.9445 0.9461 0.9461 0 0 0.9457 0.9461 0.9461 0.9461 0 0.9460 0.9461 0.9461 0.9461 0.9461

其中T为Romberg积分表,由输出结果可知计算结果为I?sinx?0xdx?0.9461.

1

例 8.5 利用Romberg积分法求I??1401?x2dx.

解:

直接利用例8.4的函数即可.

>>syms x;

>>f=sym('4/(1+x^2)') f =

4/(1+x^2)

>>[I,T]=romberg(f,0,1,5,1E-6) I = 3.1416 T =

3.0000 0 0 0 0 0 3.1000 3.1333 0 0 0 0 3.1312 3.1416 3.1421 0 0 0 3.1390 3.1416 3.1416 3.1416 0 0 3.1409 3.1416 3.1416 3.1416 3.1416 0 3.1414 3.1416 3.1416 3.1416 3.1416 3.1416

故结果为I??1401?x2dx?3.1416.

例 8.6 对积分

?1?1f(x)dx?A0f(x0)?A1f(x1)

构造其Gauss型求积公式. 解:

取f(x)?1,x,x2,x3,得到方程组:

??A0?A1?2??A0x0?A1x1?0??A0x20?A1x21?2 ?3??A330x0?A1x1?0

这是一个抽象方程组,可以利用Matlab的符号法来解之,该函数名为solve():

%直接输入解之:

>>x=solve('A0+A1=2','A0*x0+A1*x1=0','A0*x0^2+A1*x1^2=2/3','A0*x0^3+A1*x1^3=0','A0,A1,x0,x1') x =

A0: [2x1 sym] A1: [2x1 sym] x0: [2x1 sym] x1: [2x1 sym] %显示结果

>>x.A0,x.A1,x.x0,x.x1 ans = 1 1 ans = 1 1 ans =

1/3*3^(1/2) -1/3*3^(1/2) ans =

-1/3*3^(1/2) 1/3*3^(1/2)

因此得到两组解为:

(A0,A1,x0,x1)?(1,1,33,?) 3333,) 33或: (A0,A1,x0,x1)?(1,1,?求积公式为

?

1?1f(x)dx?f(33)?f(?) 33例 8.13 用中心差商公式计算f(x)?解：

应用中心差商公式可得到：

x在x?2处的一阶导数.

2?f'(2)?编制如下程序计算：

hh?2?22 h

%ex8_13.m clc;

syms x;%符号变量

f=sym('sqrt(x)')%符号函数 x0=2;%求值点

h=2;%第一次初始步长

N=10;%进行10次迭代计算 dF=diff(f,x,1)%求导

dFx=subs(dF,x0)%精确值 df=zeros(N,6); for i=1:N

df(i,1)=h;

temp=(subs(f,x,x0+h/2)-subs(f,x,x0-h/2))/h;%按中心差商公式计算 df(i,2)=double(temp); h=h/10;%步长缩小10倍 end

temp=dFx*ones(N,1)-df(:,2); df(:,3)=temp;

h=1;%第二次取初始步长为1 for i=1:N

df(i,4)=h;

temp=(subs(f,x,x0+h/2)-subs(f,x,x0-h/2))/h; df(i,5)=double(temp); h=h/10; end

temp=dFx*ones(N,1)-df(:,5); df(:,6)=temp;

disp(' 步长h 近似值 误差 步长h 近似值 误差 '); disp(df);

运行结果如下：

dFx =

0.353553390593274

步长h 近似值 误差 步长h 近似值 误差 Columns 1 through 4

2.000000000000000 0.366025403784439 -0.012472013191165 1.000000000000000 0.200000000000000 0.353663997049609 -0.000110606456335 0.100000000000000 0.020000000000000 0.353554495459696 -0.000001104866422 0.010000000000000 0.002000000000000 0.353553401641782 -0.000000011048508 0.001000000000000 0.000200000000000 0.353553390703976 -0.000000000110702 0.000100000000000

0.000020000000000 0.353553390597394 -0.000000000004120 0.000010000000000 0.000002000000000 0.353553390564087 0.000000000029186 0.000001000000000 0.000000200000000 0.353553389897954 0.000000000695320 0.000000100000000 0.000000020000000 0.353553397669515 -0.000000007076241 0.000000010000000 0.000000002000000 0.353553408771745 -0.000000018178471 0.000000001000000

Columns 5 through 6

0.356393958692601 -0.002840568099327 0.353581019507412 -0.000027628914138 0.353553666807604 -0.000000276214331 0.353553393355410 -0.000000002762136 0.353553390621819 -0.000000000028545 0.353553390586292 0.000000000006982 0.353553390564087 0.000000000029186 0.353553391008177 -0.000000000414903 0.353553386567285 0.000000004025989 0.353553408771745 -0.000000018178471

8.2 Matlab数值积分函数介绍 1. 求积分的符号运算int

I=int(fun,v,a,b)

fun为被积函数和符号表达式,可以为函数向量或函数矩阵,v为积分变量,积分函数中只有一个变量时可省略,a和b为积分上下限,即在[a,b]上计算积分,若省略则返回积分函数的一个原函数.

I为积分结果,若输出结果为符号形式的积分值时,可配合vpa(I,a)加以精度控制,其中a为精度的位数. 2. 梯形积分

I=trapz(x,y)

x和y是同维向量或矩阵,返回积分结果. 3. 复化Simpson公式数值积分quad

quad(fun,a,b,tol) quad采用自适应步长的Simpson求积法,是低阶法数值积分的函数.它自动变换、选择步长，以满足精度要求，实现变步长复合抛物形积分计算.其中，fun为积分函数,可用字符表达式,内联函数或M文件的函数形式标识,[a,b]为积分区间,tol 为积分绝对误差,默认为1E-6. 4. quadl函数

quadl(fun,a,b,tol)

采用自适应递推步长复合Lobatto 数值积分法计算积分. 5. dblquad函数

I=dblquad(fun,a,b,c,d,tol)

在矩形区域上求二重积分.fun为二元函数,可用inline定义或写成M文件函数的形式,[a,b]为变量x的上下限,[c,d]是变量y的上下限,tol 为精度要求,缺省值为1E-6. 6. triplequad函数

I=triplequad(fun,a,b,c,d,e,f,tol)

用于三重积分,在立体区域上计算三重积分结果.fun为三元函数,可用内联函数定义或写成M文件函数形式,[a,b]是变量x的积分上下限,[c,d]是变量y的积分上下限,[e,f]为变量z的积分上下限,tol为积分精度, 缺省值为1E-6. 7. 复化的8阶Newton-Cotes公式quad8

I=quad8(fun,a,b,tol) 参数意义同上. 8. 数值差分diff

D=diff(X)

这是一个非常粗略的微分函数命令diff,它计算数组间的差分.

本文来源：https://www.bwwdw.com/article/5klo.html