计算方法上机报告

计算方法上机报告

硕5028

李兰鑫

3115035005

1. 计算以下和式：0142118184858616

n

n S n n n n ∞

=??=--- ?++++??

∑

，要求： (1)若保留11个有效数字，给出计算结果，并评价计算的算法； (2)若要保留30个有效数字，则又将如何进行计算。 1.1算法思想

由题可得S<(1/16^n)*(4/(8*n+1))

保留11个有效数字，即误差ε<0.5*10^（1-11） 保留30个有效数字，即误差ε<0.5*10^（1-30） 1.2源程序

%有效数字为11的情况

digita=11;（第二问中改为30即可） for i=1:1:1000

S0=(1/(16^(i-1)))*(4/(8*i-7)); S1=(1/16^i)*(4/(8*i+1)); if ((S0-S1)<0.5*10^(1-digita)) n=i; break ; end end S=0;

digits(digita);

for i=0:1:n

S=S+vpa(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); end

1.3运行结果

2.某通信公司在一次施工中，需要在水面宽度为20

米的河沟底部沿

直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测

(1)请用合适的曲线拟合所测数据点；

(2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图；2.1算法思想

由题可知应使用多项式插值，考虑到本例中的拟合数据点较多，使用多项式插值时会因龙格现象造成误差较大，故决定采用三次样条插值。

2.2源程序

%------三次样条插值多项式------

clc;

clear;

%-----插值函数和插值点-------

n=20;%插值点个数11

i=0:n;%子元素从0到n

x=10*i;

f=[9.01 8.96 7.96 7.97 8.02 9.05 10.13 11.18 12.26 13.28 13.32 12.61 11.29 10.22

9.15 7.90 7.95 8.86 9.81 10.80 10.93];

%--输出插值点---

fprintf('当n=%d时\n',n)

fprintf('插值点x坐标为：\n')

disp(x);

fprintf('插值点y坐标f（x）为：\n')

disp(f);

%----计算差商表，并将有用的差商进行选择和保存-----

y=f;%存放差商

d=zeros(1,n+1);%存放三弯矩方程组的右端项

for i=1:n%计算第i阶差商

for j=n+1:-1:i+1

y(j)=(y(j)-y(j-1))/(x(j)-x(j-i));

end

if i==2%将第2阶差商赋值给d

d(2:n)=6.*y(3:n+1);

end

if i==3%三次样条插值的第三类边界条件

d(1)=-12*(x(2)-x(1))*y(i+1);

d(n+1)=12*(x(n+1)-x(n))*y(n+1); end

end

%-------用追赶法求三弯矩方程组------

a=ones(1,n+1);%存储下三角的非零元素b=2*ones(1,n+1);%存储对角元素

c=ones(1,n+1);%存储上三角的非零元素for i=2:n

a(i)=(x(i)-x(i-1))/(x(i+1)-x(i-1));

c(i)=1-a(i);

end

a(n+1)=-2;

c(1)=-2;

%-------计算矩阵L和U-----------

%初始化L和U

l=zeros(1,n+1);

u=zeros(1,n+1);

r=zeros(1,n+1);

u(1)=b(1);

r(1)=c(1);

for k=2:n+1

l(k)=a(k)/u(k-1);

u(k)=b(k)-l(k)*r(k-1);

r(k)=c(k);

end

%--------使用追赶法进行求解----------

M=zeros(1,n+1);%方程组的解

M(1)=d(1);

for k=2:n+1

M(k)=d(k)-l(k)*M(k-1);

end

M(n+1)=M(n+1)/u(n+1);

for k=n:-1:1

M(k)=(M(k)-r(k)*M(k+1))/u(k);

end

fprintf('三弯矩方程组的解M为：\n')

disp(M);

%------计算f（x）、N(x)和S(x)-------

Sy1=zeros(1,200);

x1=zeros(1,200);

for i=1:200

x1(i)=i-1;

end

for i=1:200

for j=1:n

if x1(i)>=x(j)&&x1(i)<=x(j+1)

temp=x(j+1)-x(j);

Sy1(i)=M(j).*(x(j+1)-x1(i)).^3/6/temp+M(j+1).*(-x(j)+x1(i)).^3/6/temp+(f(j)-M(j)*te mp^2/6).*(x(j+1)-x1(i))/temp+(f(j+1)-M(j+1)*temp^2/6).*(-x(j)+x1(i))/temp;

end

end

end

fprintf('由三次样条插值方法计算的各年的产量为：\n')

disp(Sy1)

plot(x,f,'*')%画出由插值点构成的曲线

hold on

plot(x1,Sy1,'r')%用红线画出三次样条插值拟合的曲线

hold on

legend('三次样条插值');

%------求解电缆长度------

Cordinate=[0:20/(length(Sy1)-1):20;Sy1];

diff_Cordinate=[diff(Cordinate(1,:));diff(Cordinate(2,:))];

disp('电缆长度为：');

anwser=sum(sqrt(diff_Cordinate(1,:).*diff_Cordinate(1,:)+diff(Cordinate(2,:)).*diff( Cordinate(2,:))))

2.3运行结果

由图可知，三次样条插值拟合数据点比较好，没有龙格现象。并得到，所需光缆长度约为26.6290米。

3.假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这

3.1算法思想

本题仍考虑使用三次样条插值来拟合数据点，使误差较小。

3.2源程序

X=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24];

Y=[15 14 14 14 14 15 16 18 20 20 23 25 28 31 34 31 29 27 25 24 22 20 18 17 16];

XX=0:0.01:24;

YY=spline(X,Y,XX);

coff=spline(X,Y);

plot(X,Y,'o',XX,YY);

aver=sum(YY)/length(YY); %计算区间[0,1]的误差

dp1=polyder(coff.coefs(1:4)); d2p=polyder(dp1);

x1=0:1;

f=inline(poly2sym(d2p));

y1=f(x1);

max1=max(y1);

min1=min(y1);

if(abs(max1)>abs(min1))

m1=abs(max1);

else

m1=abs(min1);

end

err1=m1/2*(0.01)^2;

3.3运行结果

上图为三次样条拟合曲线，可知拟合效果较好。并可得到，平均气温为21.4339度，误差为0.000076397。

4.设计算法，求出非线性方程52

x x

-+=的所有实根，并使误差

645200

不超过4

10-。

4.1算法思想

本题为非线性方程求解，综合考虑采用割线法求解，避免了Newton 法需要每步求解导函数的问题。

4.2源程序

%主程序

function [xx,i] = FXF(x1,x2)

e=10^(-4); %误差要求

temp=0;

i=0;

while abs(x2-x1)>=e %判断误差

temp=x2;

x2=Fx(x1,x2);

x1=temp;

i=i+1;

end

xx=x2;

end

%割线法

function [f]=Fx(x1,x2)

Fx1=6*x1^5-45*x1^2+20; Fx2=6*x2^5-45*x2^2+20;

f=x2-(x1-x2)*Fx2/Fx1/(1-Fx2/Fx1); %????·¨μü′ú1?ê? end

4．3运行结果

由题可知，一次导数有两个零点，而原函数在正负无穷大分别取到正无穷和负无穷大，因此本题共有三个实根。分别选取初值进行逼近，可以看出迭代次数相对较小。同时迭代次数还与所选初值有关，初值离真实值越远，迭代次数越多。 5.线性方程组求解。

（1）编写程序实现大规模方程组的高斯消去法程序，并对所附的方程组进行求解。所附方程组的类型为对角占优的带状方程组。 （2）针对本专业中所碰到的实际问题，提炼一个使用方程组进行求解的例子，并对求解过程进行分析、求解。 5.1算法思想

列主元高斯消去法的中心是“降维”，主要步骤包括消去和回代。针对压缩方程组，在存放时只依次存入带状部分的非零值，以提高运算速度和减少内存占用量。

(1)

(1)

()(1)(1)

()(1)(1) ,1,..., ,1,2,...,k ik

ik k kk

k k k ij ij ik kj k k k i i ik k a l a a a l a j k n

l i k k n

βββ------==-=+=-=++高斯消去法计算关系式

5.2源程序

%%读取系数矩阵

[f,p]=uigetfile('*.dat','选择数据文件'); %取数据文件 num=5; %输入系数矩阵文件头的个数 name=strcat(p,f); file=fopen(name,'r');

head=fread(file,num,'uint'); %读取二进制头文件 id=dec2hex(head(1)); fprintf('文件标识符为'); id

ver=dec2hex(head(2)); fprintf('文件版本号为'); ver

n=head(3);

fprintf('矩阵A 的阶数'); n

q=head(4);

fprintf('矩阵A 的上带宽'); q

p=head(5);

fprintf('矩阵A 的下带宽'); p

dist=4*num;

fseek(file,dist,'bof'); |

[A,count]=fread(file,inf,'float'); fclose(file);

%%对非压缩带状矩阵进行求解 if ver=='102', a=zeros(n,n); for i=1:n, for j=1:n,

a(i,j)=A((i-1)*n+j); end end

b=zeros(n,1); for i=1:n,

b(i)=A(n*n+i);

(1)(1)

(1)(1)1 ,1,2,...,1n n n n nn n

k k k k

kj j j k x a x a x k n n ββ----=+?=??

??=-=--??

∑回代过程

end

for k=1:n-1,¨

m=k;

for i=k+1:n,

if abs(a(m,k))

m=i;

end

end

if a(m,k)==0

disp('错误！')

return

end

for j=1:n,

t=a(k,j);

a(k,j)=a(m,j);

a(m,j)=t;

t=b(k);

b(k)=b(m);

b(m)=t;

end

for i=k+1:n,

a(i,k)=a(i,k)/a(k,k);

for j=k+1:n

a(i,j)=a(i,j)-a(i,k)*a(k,j);

end

b(i)=b(i)-a(i,k)*b(k);

end

end

x=zeros(n,1);

x(n)=b(n)/a(n,n);

for k=n-1:-1:1,

x(k)=(b(k)-sum(a(k,k+1:n)*x(k+1:n)))/a(k,k); end

end

%%对压缩带状矩阵进行求解

if ver=='202', ·¨

m=p+q+1;

a=zeros(n,m);

for i=1:1:n

for j=1:1:m

a(i,j)=A((i-1)*m+j);

end

end

b=zeros(n,1);

for i=1:1:n

b(i)=A(n*m+i);

end

for k=1:1:(n-1)

if a(k,(p+1))==0

disp('错误！');

break;

end

st1=n;

if (k+p)

st1=k+p;

end

for i=(k+1):1:st1

a(i,(k+p-i+1))=a(i,(k+p-i+1))/a(k,(p+1));

for j=(k+1):1:(k+q)

a(i,j+p-i+1)=a(i,j+p-i+1)-a(i,k+p-i+1)*a(k,j+p-k+1);

end

b(i)=b(i)-a(i,k+p-i+1)*b(k);

end

end

x=zeros(n,1);

x(n)=b(n)/a(n,p+1);

sum=0;

for k=(n-1):-1:1

sum=b(k);

st2=n;

if (k+q)

st2=k+q;

end

for j=(k+1):1:st2

sum=sum-a(k,j+p-k+1)*x(j);

end

x(k)=sum/a(k,p+1);

sum=0;

end

end

disp('方程组的的解为：')

disp(x)

5.3运行结果

打开dat51、dat52后的结果如下图所示，分别显示了矩阵的阶数、上下带宽以及得到的解，与题目中所给出的解相同，故可以证明程序

的正确性。

分别打开dat53、dat54文件，得到x解如下图所示（由于矩阵阶数过高，只截取部分结果）。