数值分析 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;

https://www.360docs.net/doc/452006522.html,

-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求积分值：

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

sym x

int(4*7782.5*sqrt(1-(972.5/7782.5)^2*sin(x).^2),0,pi/2)

I=vpa(int

(4*7782.5*sqrt(1-(972.5/7782.5)^2*sin(x).^2),0,pi/2))

回车后输出积分结果为：

I =48707.60111

下面编写M文件用复合的辛普森公式计算周长：

function Simpson4

e=7782.5 ; %因为下面要用到a，为了避免变量冲突，所以这里将

c=972.5; %a改成e了。

x=0:pi/16:pi/2;

f=4*e*sqrt(1-(c/e)^2*(sin(x)).^2);

a=0,b=pi/2;

h=(b-a)/4;

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

for k=1:1:4

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

end

for k=2:1:4

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

end

s=(h/6)*s

运行上述程序得出周长为：4.870743851190015e+004km.

用梯形求积分函数trapz求积分

在命令行输入下列程序：

d=pi/16;

x=0:d:pi/2;

T=d*trapz(4*7782.5*sqrt(1-(972.5/7782.5)^2*sin(x).^2))

回车后，可见 T= 4.870743851190014e+004km

用复化的梯形公式求解该题：

Matlab程序如下：

function t=tixing()%梯形法求积分

clear;

clc;

options={'积分下限a','积分上限b' ,'插入点相关的值n'};

topic='seting';

lines=1;

def={'0','1','1000'};

h=inputdlg(options,topic,lines,def);

a=eval(h{1});%积分下限

b=eval(h{2});%积分上限

n=eval(h{3});%子区间个数

%********************************************

f='func';%用f来调用被积函数func

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

t=0;

for k=1:1:n-1

x=a+h*(k);

t=t + feval(f,x);

end

format long

t=h*( feval(f,a) +feval(f,b)+2*t)/2;%s是梯形规则的总计

end

%定义被积函数func

function y=func(x)

y=4*7782.5*sqrt(1-(972.5/7782.5)^2*sin(x).^2);

end

运行上面的程序得周长为：4.870743851190014e+004。

由上述计算结果可知编程计算与用积分函数计算的结果是一样的。与精确值之间的误差也在允许的范围内。

终上所述，卫星轨道的周长约为48707km.