数值分析 matlab 实验4

(1) 解题过程如下：

（1）MATLAB中创建复化梯形公式和复化辛普森公式的 M 文件：1）复化梯形公式文件：

function s=T_fuhua(f,a,b,n)

h=(b-a)/n;

s=0;

for k=1:(n-1)

x=a+h*k;

s=s+feval(f,x);

end

s=h*(feval(f,a)+feval(f,b))/2+h*s;

2）复化辛普森公式文件：

function s=S_fuhua(f,a,b,n)

h=0;

h=(b-a)./(2*n);

s1=0;

-5-

s2=0;

for k=1:n-1

x=a+h*2*k;

s1=s1+feval(f,x);

end

for k=1:n

x=a+h*(2*k-1);

s2=s2+feval(f,x);

end

s=h*(feval(f,a)+feval(f,b)+s1*2+s2*4)/3;

在MATLAB中输入：

f=inline('x/(4+x^2)');a=0;b=1;

%inline 构造内联函数对象

for n=2:10

s(n-1)=T_fuhua(f,a,b,n);s(n-1)=vpa(s(n-1),10);

%调用复化梯形公式，生成任意精度的数值

end

exact=int('x/(4+x^2)',0,1);exact=vpa(exact,10)

%求出积分的精确值

输出结果：exact =

.1115717755

s =

Columns 1 through 6

0.1088 0.1104 0.1109 0.1111 0.1113 0.1114

Columns 7 through 9

0.1114 0.1114 0.1115

在MATLAB中输入以下函数用以画出计算误差与 n 之间的曲线：

r=abs(exact-s);

n=2:10;

plot(double(n),double(r(n-1)))

得到结果如图所示：

（2）在 MATLAB中输入以下程序代码：

f=inline('x/(4+x^2)');a=0;b=1;n=9;

%inline 构造内联函数对象

t=T_fuhua(f,a,b,n);t=vpa(t,10)

s=S_fuhua(f,a,b,n);s=vpa(s,10)

%调用复化梯形和复化辛普森公式，生成任意精度的数值

exact=int('x/(4+x^2)',0,1);exact=vpa(exact,10)

%求出积分的精确值

计算结果：t =

.1114379370

s =

.1115717991

exact =

.1115717755

E1=|t-exact|=0.0001338385

E2=|s-exact|=0.0000000236

所以,两种方法计算所得的绝对误差:E1>E2

（1）中的两个结果 s 与t，两个函数的计算量基本相同，但是精度却有很大

差别：与精确值exact =.1115717755比较，复化梯形公式的结果t =.1114379370 只有三位有

效数字，而复化辛普森公式的结果 s =.1115717991 却有七位有效数字。

由此例子可以知道，复化辛普森公式比复化梯形公式的代数精度高，更适合于科学计

算与应用。

(2)用求定积分的函数int求积分的值：

在命令窗口中键入如下的命令：

sym x

int(sqrt(x),1,9)

I=vpa(int (sqrt(x),1,9))

回车后输出积分结果为：

I =17.3333

编写求复合辛普森公式的m文件如下：

function Simpson2

x=1:1/2:9;

f= sqrt(x);

a=1,b=9;

%f(1)=1;

h=(b-a)/4;

s=f(1)+f(17);

for k=1:1:8

s=s+4*f(2*k);

end

for k=2:1:8

s=s+2*f(2*k-1);

end

s=(h/6)*s

运行上述程序：s=17.3332

(2)用求定积分的函数int求积分的值：

在命令窗口中键入如下的命令：

sym x

int(sqrt(4-sin(x).^2),0,pi/6)

I=vpa(int (sqrt(4-sin(x).^2),0,pi/6))

回车后输出积分结果为：

I =1.03576

编写求复合辛普森公式的m文件如下：

function Simpson3

x=0:pi/96:pi/6;

f=sqrt(4-sin(x).^2);

a=0,b=pi/6;

%f(1)=1;

h=(b-a)/6;

s=f(1)+f(17);

for k=1:1:8

s=s+4*f(2*k);

end

for k=2:1:8

s=s+2*f(2*k-1);

end

s=(h/6)*s

运行上述程序：s= 1.0358 误差为0.00004

分析：本题是求解定积分的问题，可以采用复化的辛普森公式求解。根据题中的条件可以计算出c,a的值。

解：用matlab求解才c,a如下：

a =(2R+H+h)/2 ,c=(H-h)/2. c=972.5, a=7782.5

用求定积分的函数int求积分值：