数值分析 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求积分值：
在命令窗口中键入如下的命令：
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.。