西安交通大学 计算方法上机报告

计算方法上机报告

姓名：

学号：

班级：机械硕4002

上课班级：02班

说明：

本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。

1. 对以下和式计算：

∑

∞

?

?? ??+-+-+-+=0681581482184161n n n n S n

，要求：

① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算；

（1） 算法思想

1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114

16818485861681

n n n

a n n n n n ε??=

---<< ?+++++??； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数)

（2）算法结构

1. ;0=s

??

?

??+-+-+-+=

681581482184161n n n n t n

； 2. for 0,1,2,,n i =??? if 10m t -≤

end;

3. for ,1,2,,0n i i i =--???

;s s t =+

（3）Matlab源程序

clear; %清除工作空间变量

clc; %清除命令窗口命令

m=input('请输入有效数字的位数m='); %输入有效数字的位数

s=0;

for n=0:50

t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6));

if t<=10^(-m) %判断通项与精度的关系break;

end

end;

fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值

for i=n-1:-1:0

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

s=s+t; %求和运算

end

s=vpa(s,m) %控制s的精度

（4）结果与分析

当保留11位有效数字时，需要将n值加到n=7，

s =3.1415926536；

当保留30位有效数字时，需要将n值加到n=22,

s =3.14159265358979323846264338328。

通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。

2.某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线

走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：

分点0 1 2 3 4 5 6

深度9.01 8.96 7.96 7.97 8.02 9.05 10.13

分点7 8 9 10 11 12 13

深度11.18 12.26 13.28 13.32 12.61 11.29 10.22

分点14 15 16 17 18 19 20

深度9.15 7.90 7.95 8.86 9.81 10.80 10.93

①请用合适的曲线拟合所测数据点；

②预测所需光缆长度的近似值，作出铺设河底光缆的曲线图；

（1）算法思想

如果使用多项式差值，则由于龙格现象，误差较大，因此，用相对较少的插值数据点作插值，可以避免大的误差，但是如果又希望将所得数据点都用上，且所用数据点越多越好，可以采用分段插值方式，即用分段多项式代替单个多项式作插值。分段多项式是由一些在相互连接的区间上的不同多项式连接而成的一条连续曲线，其中三次样条插值方法是一种具有较好“光滑性”的分段插值方法。

在本题中，假设所铺设的光缆足够柔软，在铺设过程中光缆触地走势光滑，紧贴地面，并且忽略水流对光缆的冲击。海底光缆线的长度预测模型如下所示，光缆从A点铺

h。

至B点，在某点处的深度为i

海底光缆线的长度预测模型

计算光缆长度时，用如下公式：

20

()L f x ds =

?

20

'20()1()f x f x dx =+? 19

1'20

()1()k k

k f x f x dx

+==+∑?

22

()()x y =?+?∑

（2）算法结构

1. For n i ,,2,1,0???=

1.1 i i M y ?

2. For 2,1=k

2.1 For k n n i ,,1, -=

2.1.1 i k i i i i M x x M M ?----)/()(1

3. 101h x x ?-

4. For 1-,,2,1n i =

4.1 11++?-i i i h x x

4.2 b a c c h h h i i i i i i ??-?+++2;1;)/(11

4.3 i i d M ?+16

5. 0000;;c M d M d n n ???λ

n n n b a b ???2;;20μ

6. 1111,γμ??d b

7. 获取M 的矩阵元素个数，存入m 8. For m k ,,3,2 =

8.1

k k k l a ?-1/μ

8.2 k k k k c l b μ??-1- 8.3 k k k k l d γγ??-1- 9. m m m M ?μγ/

10. For 1,,2,1 --=m m k

10.1 k k k k k M M c ??-+μγ/)(1 11. 获取x 的元素个数存入s 12. k ?1

13. For 1,,2,1-=s i

13.1 if i x x ≤~ then k i ?；break

else k i ?+1

14. x

x x x x x h x x k k k k ?~;~;

11?-?-?--- y h x h M y x h M y x M x M k k k k k k ~/]?)6

()6(6?6[2

211331

?-+-++---

（3）Matlab 源程序

clear; clc;

x=0:1:20; %产生从0到20含21个等分点的数组

X=0:0.2:20;

y=[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]; %等分点位置的深度数据

n=length(x); %等分点的数目

N=length(X);

%% 求三次样条插值函数s(x)

M=y;

for k=2:3; %计算二阶差商并存放在M中for i=n:-1:k;

M(i)=(M(i)-M(i-1))/(x(i)-x(i-k+1));

end

end

h(1)=x(2)-x(1); %计算三对角阵系数a,b,c及右端向量d for i=2:n-1;

h(i)=x(i+1)-x(i);

c(i)=h(i)/(h(i)+h(i-1));

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

b(i)=2;

d(i)=6*M(i+1);

end

M(1)=0; %选择自然边界条件

M(n)=0;

b(1)=2;

b(n)=2;

c(1)=0;

a(n)=0;

d(1)=0;

d(n)=0;

u(1)=b(1); %对三对角阵进行LU分解y1(1)=d(1);

for k=2:n;

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

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

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

end

M(n)=y1(n)/u(n); %追赶法求解样条参数M(i) for k=n-1:-1:1;

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

end

s=zeros(1,N);

for m=1:N;

k=1;

for i=2:n-1

if X(m)<=x(i);

k=i-1;

break;

else

k=i;

end

end

H=x(k+1)-x(k); %在各区间用三次样条插值函数计算X点处的值x1=x(k+1)-X(m);

x2=X(m)-x(k);

s(m)=(M(k)*(x1^3)/6+M(k+1)*(x2^3)/6+(y(k)-(M(k)*(H^2)/6))*x1+(y(k+1)-(M(k+1)*(H^2)/ 6))*x2)/H;

end

%% 计算所需光缆长度

L=0; %计算所需光缆长度

for i=2:N

L=L+sqrt((X(i)-X(i-1))^2+(s(i)-s(i-1))^2);

end

disp('所需光缆长度为L=');

disp(L);

figure

plot(x,y,'*',X,s,'-') %绘制铺设河底光缆的曲线图xlabel('位置','fontsize',16); %标注坐标轴含义

ylabel('深度/m','fontsize',16);

title('铺设河底光缆的曲线图','fontsize',16);

grid;

（4）结果与分析

铺设海底光缆的曲线图如下图所示：

仿真结果表明，运用分段三次样条插值所得的拟合曲线能较准确地反映铺设光缆的走势图，计算出所需光缆的长度为L=26.4844m。

3. 假定某天的气温变化记录如下表所示，试用数据拟合的方法找出这一天的气温变化的规律；试计算这一天的平均气温，并试估计误差。

时刻0 1 2 3 4 5 6 7 8 9 10 11 12

平均气温15 14 14 14 14 15 16 18 20 20 23 25 28

时刻13 14 15 16 17 18 19 20 21 22 23 24

平均气温31 34 31 29 27 25 24 22 20 18 17 16

（1）算法思想

在本题中，数据点的数目较多。当数据点的数目很多时，用“多项式插值”方法做数据近似要用较高次的多项式，这不仅给计算带来困难，更主要的缺点是误差很大。用“插值样条函数”做数据近似，虽然有很好的数值性质，且计算量也不大，但存放参数Mi 的量很大，且没有一个统一的数学公式来表示，也带来了一些不便。另一方面，在有的实际问题中，用插值方法并不合适。当数据点的数目很大时，要求)(x p 通过所有数据点，可能会失去原数据所表示的规律。如果数据点是由测量而来的，必然带有误差，插值法要求准确通过这些不准确的数据点是不合适的。在这种情况下，不用插值标准而用其他近似标准更加合理。通常情况下，是选取i α使2E 最小，这就是最小二乘近似问题。

在本题中，采用“最小二乘法”找出这一天的气温变化的规律，使用二次函数、三次函数、四次函数以及指数型函数2

)(c t b ae C --=，计算相应的系数，估算误差，并作图比较各种函数之间的区别。

（2）算法结构

本算法用正交化方法求数据的最小二乘近似。假定数据以用来生成了G ，并将

y 作为其最后一列（第1+n 列）存放。结果在a 中，ρ是误差22E 。

I 、使用二次函数、三次函数、四次函数拟合时

1. 将“时刻值”存入

i x ，数据点的个数存入m

2. 输入拟合多项式函数)(x p 的最高项次数1-=n h ，则拟合多项式函数为

)()()()(2211x g x g x g x p n n ααα+???++= , 根据给定数据点确定G For 1,,2,1,0-???=n j

For m i ,,2,1???= 2.1 1,+?j i j

i g x 2.2 1,+?n i i g y

3. For n k ,,2,1???=

3.1 [形成矩阵k Q ]

3.1.1 σ?-∑=m

k i ik

kk g g 2

/12))(sgn( 3.1.2 k kk g ωσ?-

3.1.3 For m k k j ,,2,1???++=

3.1.3.1 j jk g ω? 3.1.4 βσω?k 3.2 [变换1-k G 到k G ] 3.2.1 kk g ?σ

For 1,,,2,1+???++=n n k k j 3.2.2 t g m

k i ij i ?∑=βω/)(

3.2.3 For m k k i ,,1,???+=

3.2.3.1 ij i ij g t g ?+ω

4. [解三角方程1h Ra =] 4.1

n nn n n a g g ?+/1.

4.2 For 1,,2,1???--=n n i 4.2.1

i

ii n

i j j ij

n i a g x g

g ?-

∑+=+/][1

1.

5. [计算误差2

2E ]

ρ

?∑+=+m

n i n i g 1

2

1,

II 、使用指数函数拟合时

现将指数函数进行变形：

将y C =,x t =代入2

)(c t b ae

C --=得:

2

)(c x b ae

y --=

对上式左右取对数得: 2

2

2ln ln bx bcx bc a y -+-= 令b bc bc a y z -==-==212

02ln ln ααα，，，

则可得多项式: 2

210x x z ααα++=

（3）Matlab 源程序

clear; %清除工作空间变量 clc; %清除命令窗口命令 x=0: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]; [~,m]=size(x); %将数据点的个数存入m T=sum(y(1:m))/m;

fprintf('一天的平均气温为T=%f\n',T); %求一天的平均气温 %% 二次、三次、四次函数的最小二乘近似

h=input('请输入拟合多项式的最高项次数='); %根据给定数据点生成矩阵G n=h+1; G=[];

for j=0:(n-1)

g=x.^j; %g(x)按列排列 G=vertcat(G,g); %g 垂直连接G end

G=G'; %转置得到矩阵G

for i=1:m %将数据y 作为G 的最后一列（n+1列）

G(i,n+1)=y(i);

end

G;

for k=1:n %形成矩阵Q(k) if G(k,k)>0;

sgn=1;

elseif G(k,k)==0;

sgn=0;

else sgn=-1;

end

sgm=-sgn*sqrt(sum(G(k:m,k).^2));

w=zeros(1,n);

w(k)=G(k,k)-sgm;

for j=k+1:m

w(j)=G(j,k);

end

bt=sgm*w(k);

G(k,k)=sgm; %变换Gk-1到Gk for j=k+1:n+1

t=sum(w(k:m)*G(k:m,j))/bt;

for i=k:m;

G(i,j)=G(i,j)+t*w(i);

end

end

end

A (n)=G(n,n+1)/G(n,n); %解三角方程求系数A

for i=n-1:-1:1

A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);

end

e=sum(G(n+1:m,n+1).^2); %计算误差e

fprintf('%d次函数的系数是：',h); %输出系数a及误差e disp(A);

fprintf('使用%d次函数拟合的误差是%f：',h,e);

t=0:0.05:24;

A=fliplr(A); %将系数数组左右翻转

Y=poly2sym(A); %将系数数组转化为多项式subs(Y,'x',t);

Y=double(ans);

figure(1)

plot(x,y,'k*',t,Y,'r-'); %绘制拟合多项式函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');

title([num2str(n-1),'次函数的最小二乘曲线']);

grid;

%% 指数函数的最小二乘近似

yy=log(y);

n=3;

G=[];

GG=[];

for j=0:(n-1)

g=x.^j; %g(x)按列排列

G=vertcat(G,g); %g垂直连接G

gg=t.^j; %g(x)按列排列

GG=vertcat(GG,gg); %g垂直连接G

end

G=G'; %转置得到矩阵G

for i=1:m %将数据y作为G的最后一列（n+1列）G(i,n+1)=yy(i);

end

G;

for k=1:n %形成矩阵Q(k)

if G(k,k)>0;

sgn=1;

elseif G(k,k)==0;

sgn=0;

else sgn=-1;

end

sgm=-sgn*sqrt(sum(G(k:m,k).^2));

w=zeros(1,n);

w(k)=G(k,k)-sgm;

for j=k+1:m

w(j)=G(j,k);

end

bt=sgm*w(k);

G(k,k)=sgm; %变换Gk-1到Gk

for j=k+1:n+1

t=sum(w(k:m)*G(k:m,j))/bt;

for i=k:m;

G(i,j)=G(i,j)+t*w(i);

end

end

end

A(n)=G(n,n+1)/G(n,n); %解三角方程求系数A for i=n-1:-1:1

A (i)=(G(i,n+1)-sum(G(i,i+1:n).*A (i+1:n)))/G(i,i);

end

b=-A(3);

c=A(2)/(2*b);

a=exp(A(1)+b*(c^2));

G(n+1:m,n+1)=exp(sum(G(n+1:m,n+1).^2));

e=sum(G(n+1:m,n+1).^2); %计算误差e

fprintf('\n指数函数的系数是：a=%f,b=%f,c=%f',a,b,c); %输出系数及误差e fprintf('\n使用指数函数拟合的误差是：%f',e);

t=0:0.05:24;

YY=a.*exp(-b.*(t-c).^2);

figure(2)

plot(x,y,'k*',t,YY,'r-'); %绘制拟合指数函数图形xlabel('时刻'); %标注坐标轴含义ylabel('平均气温');

title(['指数函数的最小二乘曲线']);

grid;

（4）结果与分析

a、二次函数：

一天的平均气温为：21.2000

2次函数的系数: 8.3063 2.6064 -0.0938

使用2次函数拟合的误差是：280.339547

二次函数的最小二乘曲线如下图所示：

b、三次函数：

一天的平均气温为：21.2000

3次函数的系数: 13.3880 -0.2273 0.2075 -0.0084 使用3次函数拟合的误差是：131.061822

三次函数的最小二乘曲线如下图所示：

c、四次函数：

一天的平均气温为：21.2000

4次函数的系数: 16.7939 -3.7050 0.8909 -0.0532 0.0009 使用4次函数拟合的误差是：59.04118

四次函数的最小二乘曲线如下图所示：

d、指数函数：

一天的平均气温为：21.2000

指数函数的系数是：a=26.160286，b=0.004442，c=14.081900

使用指数函数拟合的误差是：57.034644

指数函数的最小二乘曲线如下图所示：

《计算方法》课内实验报告

《计算方法》实验报告 姓名： 班级： 学号： 实验日期: 2011年10月26日

一、实验题目： 数值积分 二、实验目的： 1．熟悉matlab 编写及运行数值计算程序的方法。 2．进一步理解数值积分的基础理论。 3．进一步掌握应用不同的数值积分方法求解给定的积分并给出数据结果及误差分析。 三、实验内容： 1．分别用复合梯形求积公式及复合辛普森求积公式计算积分xdx x ln 10 ? , 要求计算精度达到410-,给出计算结果并比较两种方法的计算节点数. 2．用龙贝格求积方法计算积分dx x x ?+3 021,使误差不超过510-. 3．用3=n 的高斯-勒让德公式计算积分?3 1 sin x e x ,给出计算结果. 4．用辛普森公式(取2==M N ) 计算二重积分.5 .00 5 .00 dydx e x y ? ? - 四、实验结果： 1.（1）复合梯形法: 将区间[a,b]划分为n 等份，分点n k n a b h kh a x k ,2,1,0,,=-=+=在每个区间[1,+k k x x ](k=0,1,2,···n-1)上采用梯形公式，则得 )()]()([2)()(1 11 1 f R x f x f h dx x f dx x f I n n k k k b a n k x x k k ++===∑?∑? -=+-=+ 故)]()(2)([21 1 b f x f a f h T n k k n ++=∑-=称为复合梯形公式 计算步长和划分的区间 Eps=1E-4 h1=sqrt(Eps/abs(-(1-0)/12*1/(2+1))) h1 =0.0600 N1=ceil(1/h1) N1 =17 用复合梯形需要计算17个结点。 复合梯形： function T=trap(f,a,b,n) h=(b-a)/n;

统计西安交大期末考试试题(含答案)

西安交大统计学考试试卷 一、单项选择题（每小题2 分，共20 分） 1.在企业统计中，下列统计标志中属于数量标志的是（C） A、文化程度 B、职业 C、月工资 D、行业 2.下列属于相对数的综合指标有（B ） A、国民收入 B、人均国民收入 C、国内生产净值 D、设备台数 3.有三个企业的年利润额分别是5000 万元、8000 万元和3900 万元，则这句话中有（B）个变量？ A、0 个 B、两个 C、1 个 D、3 个 4.下列变量中属于连续型变量的是（A ） A、身高 B、产品件数 C、企业人数 D、产品品种 5.下列各项中，属于时点指标的有（A ） A、库存额 B、总收入 C、平均收入 D、人均收入 6.典型调查是（B ）确定调查单位的 A、随机 B、主观 C、随意 D 盲目 7.总体标准差未知时总体均值的假设检验要用到（A ）： A、Z 统计量 B、t 统计量 C、统计量 D、X 统计量 8.把样本总体中全部单位数的集合称为（A ） A、样本 B、小总体 C、样本容量 D、总体容量 9.概率的取值范围是p（D ） A、大于1 B、大于－1 C、小于1 D、在0 与1 之间 10.算术平均数的离差之和等于（A ） A、零 B、1 C、－1 D、2 二、多项选择题（每小题2 分，共10 分。每题全部答对才给分，否则不计分） 1.数据的计量尺度包括（ABCD ）： A、定类尺度 B、定序尺度 C、定距尺度 D、定比尺度 E、测量尺度 2.下列属于连续型变量的有（BE ）： A、工人人数 B、商品销售额 C、商品库存额 D、商品库存量 E、总产值 3.测量变量离中趋势的指标有（ABE ） A、极差 B、平均差 C、几何平均数 D、众数 E、标准差 4.在工业企业的设备调查中（BDE ） A、工业企业是调查对象 B、工业企业的所有设备是调查对象 C、每台设备是 填报单位D、每台设备是调查单位E、每个工业企业是填报单位 5.下列平均数中，容易受数列中极端值影响的平均数有（ABC ） A、算术平均数 B、调和平均数 C、几何平均数 D、中位数 E、众数 三、判断题（在正确答案后写“对”，在错误答案后写“错”。每小题1 分，共10 分） 1、“性别”是品质标志。（对）

计算方法上机实验报告

. / 《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：XXX 二〇一八年五月二十五日

前言 通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。 以下为本次上机实验报告，按照实验内容共分为六部分。 实验一： 一、实验名称及题目： Newton 迭代法 例2.7(P38):应用Newton 迭代法求在附近的数 值解，并使其满足. 二、解题思路： 设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交

点的横坐标) (') (0001x f x f x x - =,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标) (') (1112x f x f x x - =称2x 为'x 的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把) (') (1n n n n x f x f x x - =+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。 三、Matlab 程序代码： function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0); y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1; while abs(x1-x0)>=tol x0=x1; y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; end x=double(x1) K 四、运行结果：

数值分析上机作业

数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：

课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m )

则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。

太原理工大学数值计算方法实验报告

本科实验报告 课程名称：计算机数值方法 实验项目：方程求根、线性方程组的直接解 法、线性方程组的迭代解法、代数插值和最 小二乘拟合多项式 实验地点：行勉楼 专业班级： ******** 学号： ********* 学生姓名： ******** 指导教师：李誌，崔冬华 2016年 4 月 8 日

y = x*x*x + 4 * x*x - 10; return y; } float Calculate(float a,float b) { c = (a + b) / 2; n++; if (GetY(c) == 0 || ((b - a) / 2) < 0.000005) { cout << c <<"为方程的解"<< endl; return 0; } if (GetY(a)*GetY(c) < 0) { return Calculate(a,c); } if (GetY(c)*GetY(b)< 0) { return Calculate(c,b); } } }; int main() { cout << "方程组为：f(x)=x^3+4x^2-10=0" << endl; float a, b; Text text; text.Getab(); a = text.a; b = text.b; text.Calculate(a, b); return 0; } 2.割线法： // 方程求根（割线法）.cpp : 定义控制台应用程序的入口点。// #include "stdafx.h" #include"iostream"

心得体会 使用不同的方法，可以不同程度的求得方程的解，通过二分法计算的程序实现更加了解二分法的特点，二分法过程简单，程序容易实现，但该方法收敛比较慢一般用于求根的初始近似值，不同的方法速度不同。面对一个复杂的问题，要学会简化处理步骤，分步骤一点一点的循序处理，只有这样，才能高效的解决一个复杂问题。

西安交通大学计算方法B上机试题

1.计算以下和式：01421181 84858616n n S n n n n ∞ =?? =--- ?++++??∑ ，要求： (1)若保留11个有效数字，给出计算结果，并评价计算的算法； (2)若要保留30个有效数字，则又将如何进行计算。 （1）题目分析 该题是对无穷级数求和，因此在使用matlab 进行累加时需要一个累加的终止条件。这里令?? ? ??+-+-+-+= 681581482184161n n n n a n n ，则 ()()1.016 1 6855844864816114851384128698161 681581482184161148113811282984161111<< ? ??? ????? ??++++++???? ????? ??++++++=??? ????? ??+-+-+-+??? ????? ??+-+-+-+=+++n n n n n n n n n n n n n n n n a a n n n n n n 故近似取其误差为1+≈k a ε，并且有m -1m -111021 21 ?=?=≈+βεk a ， （2）算法依据 使用matlab 编程时用digits 函数和vpa 函数来控制位数。 （3）Matlab 运行程序 %%保留11位有效数字 k1=11; s1=0;%用于存储这一步计算值 for n=0:50 a=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); n1=n-1; if a<=0.5*10^(1-k1) break end end; for i=0:1:n1 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s1=s1+t; end s11=vpa(s1,k1); disp('保留11位有效数字的结果为：');disp(s11); disp('此时n 值为：');disp(n1); %%保留30位有效数字 clear all; k2=30;

西安交通大学计算方法B上机报告

计算方法上机报告

姓名： 学号： 班级：能动上课班级：

题目及求解： 一、对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=0681581482184161n n n n S n ，要求： ① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算； 1 算法思想 （1）根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114 16818485861681 n n n a n n n n n ε??= ---<< ?+++++??； （2）为了保证计算结果的准确性，写程序时，从后向前计算； （3）使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数) 2 算法结构 ;0=s ?? ? ??+-+-+-+= 681581482184161n n n n t n ； for 0,1,2,,n i =??? if 10m t -≤ end; for ,1,2,,0n i i i =--??? ;s s t =+ 3 Matlab 源程序 clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0;

for n=0:50 t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); if t<=10^(-m) %判断通项与精度的关系break; end end; fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值 for i=n-1:-1:0 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s=s+t; %求和运算 end s=vpa(s,m) %控制s的精度 4 结果与分析 若保留11位有效数字，则n=7，此时求解得： s =3.1415926536； 若保留30位有效数字时，则n=22, 此时求解得： s =3.8。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。 二、某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示：

数值分析上机题目详解

第一章 一、题目 设∑ =-= N N j S 2 j 2 1 1，其精确值为)11 123(21+--N N 。 1) 编制按从大到小的顺序1 1 13112122 2-+??+-+-=N S N ，计算S N 的通用程序。 2) 编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 3) 按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） 4) 通过本次上机题，你明白了什么？ 二、通用程序 N=input('Please Input an N (N>1):'); AccurateValue=single((0-1/(N+1)-1/N+3/2)/2); Sn1=single(0); for a=2:N; Sn1=Sn1+1/(a^2-1); end Sn2=single(0); for a=2:N; Sn2=Sn2+1/((N-a+2)^2-1); end fprintf('The value of Sn (N=%d)\n',N); fprintf('Accurate Calculation %f\n',AccurateValue); fprintf('Caculate from large to small %f\n',Sn1); fprintf('Caculate from small to large %f\n',Sn2); disp('____________________________________________________')

三、结果 从结果可以看出有效位数是6位。 感想：可以得出，算法对误差的传播有一定的影响，在计算时选一种好的算法可以使结果更为精确。从以上的结果可以看到从大到小的顺序导致大数吃小数的现象，容易产生较大的误差，求和运算从小数到大数所得到的结果才比较准确。

c 计算器实验报告

简单计算器 姓名: 周吉祥 实验目的：模仿日常生活中所用的计算器，自行设计一个简单的计算器程序，实现简单的计算功能。 实验内容： (1)体系设计： 程序是一个简单的计算器，能正确输入数据，能实现加、减、乘、除等算术运算，运算结果能正确显示，可以清楚数据等。 (2)设计思路： 1)先在Visual C++ 6.0中建立一个MFC工程文件，名为 calculator. 2)在对话框中添加适当的编辑框、按钮、静态文件、复选框和单 选框 3)设计按钮，并修改其相应的ID与Caption. 4)选择和设置各控件的单击鼠标事件。 5)为编辑框添加double类型的关联变量m_edit1. 6)在calculatorDlg.h中添加math.h头文件，然后添加public成 员。 7)打开calculatorDlg.cpp文件，在构造函数中，进行成员初始 化和完善各控件的响应函数代码。 (3)程序清单：

●添加的public成员： double tempvalue; //存储中间变量 double result; //存储显示结果的值 int sort; //判断后面是何种运算：1.加法2.减法3. 乘法 4.除法 int append; //判断后面是否添加数字 ●成员初始化： CCalculatorDlg::CCalculatorDlg(CWnd* pParent /*=NULL*/) : CDialog(CCalculatorDlg::IDD, pParent) { //{{AFX_DATA_INIT(CCalculatorDlg) m_edit1 = 0.0; //}}AFX_DATA_INIT // Note that LoadIcon does not require a subsequent DestroyIcon in Win32 m_hIcon = AfxGetApp()->LoadIcon(IDR_MAINFRAME); tempvalue=0; result=0; sort=0; append=0; }

计算方法上机实验报告

《计算方法》上机实验报告 班级：XXXXXX 小组成员：XXXXXXX XXXXXXX XXXXXXX XXXXXXX 任课教师：XXX 二〇一八年五月二十五日

前言 通过进行多次的上机实验，我们结合课本上的内容以及老师对我们的指导，能够较为熟练地掌握Newton 迭代法、Jacobi 迭代法、Gauss-Seidel 迭代法、Newton 插值法、Lagrange 插值法和Gauss 求积公式等六种算法的原理和使用方法，并参考课本例题进行了MATLAB 程序的编写。 以下为本次上机实验报告，按照实验内容共分为六部分。 实验一： 一、实验名称及题目： Newton 迭代法 例2.7(P38):应用Newton 迭代法求 在 附近的数值解 ，并使其满足 . 二、解题思路： 设'x 是0)(=x f 的根,选取0x 作为'x 初始近似值,过点())(,00x f x 做曲线)(x f y =的切线L ,L 的方程为))((')(000x x x f x f y -+=,求出L 与x 轴交点的横坐标) (') (0001x f x f x x - =,称1x 为'x 的一次近似值,过点))(,(11x f x 做曲线)(x f y =的切线,求该切线与x 轴的横坐标) (') (1112x f x f x x - =称2x 为'x

的二次近似值,重复以上过程,得'x 的近似值序列{}n x ,把 ) (') (1n n n n x f x f x x - =+称为'x 的1+n 次近似值，这种求解方法就是牛顿迭代法。 三、Matlab 程序代码： function newton_iteration(x0,tol) syms z %定义自变量 format long %定义精度 f=z*z*z-z-1; f1=diff(f);%求导 y=subs(f,z,x0); y1=subs(f1,z,x0);%向函数中代值 x1=x0-y/y1; k=1; while abs(x1-x0)>=tol x0=x1; y=subs(f,z,x0); y1=subs(f1,z,x0); x1=x0-y/y1;k=k+1; end x=double(x1) K 四、运行结果： 实验二：

数值计算方法I上机实验考试题

数值计算方法I 上机实验考试题（两题任选一题） １．小型火箭初始质量为900千克，其中包括600千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生30000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数为0.4（千克/米）．重力加速度取9.8米/秒2. A. 建立火箭升空过程的数学模型（微分方程）； B. 求引擎关闭瞬间火箭的高度、速度、加速度，及火箭到达最高点的时间和高度． ２．小型火箭初始质量为1200千克，其中包括900千克燃料。火箭竖直向上发射时燃料以15千克/秒的速率燃烧掉，由此产生40000牛顿的恒定推力．当燃料用尽时引擎关闭。设火箭上升的整个过程中，空气阻力与速度平方成正比，比例系数记作k ，火箭升空过程的数学模型为 0)0(,0,01222==≤≤-+?? ? ??-==t dt dx x t t mg T dt dx k dt x d m 其中)(t x 为火箭在时刻t 的高度，m =1200-15t 为火箭在时刻t 的质量，T （=30000牛顿）为推力，g (=9.8米/秒2)为重力加速度, t 1 (=900/15=60秒)为引擎关闭时刻． 今测得一组数据如下（t ~时间（秒），x ~高度（米），v ~速度（米/秒））： 现有两种估计比例系数k 的方法： 1．用每一个数据(t,x,v )计算一个k 的估计值（共11个），再用它们来估计k 。 2．用这组数据拟合一个k ． 请你分别用这两种方法给出k 的估计值，对方法进行评价，并且回答，能否认为空气阻力系数k=0.5（说明理由）．

计算方法实验报告格式

计算方法实验报告格式 小组名称: 组长姓名(班号)： 小组成员姓名(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一个完整的实验，应包括数据准备、理论基础、实验内容及方法，最终对实验结果进行分析，以达到对理论知识的感性认识，进一步加深对相关算法的理解，数值实验以实验报告形式完成，实验报告格式如下： 一、实验名称 实验者可根据报告形式需要适当写出． 二、实验目的及要求 首先要求做实验者明确，为什么要做某个实验，实验目的是什么，做完该实验应达到什么结果，在实验过程中的注意事项，实验方法对结果的影响也可以以实验目的的形式列出． 三、算法描述（实验原理与基础理论） 数值实验本身就是为了加深对基础理论及方法的理解而设置的，所以要求将实验涉及到的理论基础，算法原理详尽列出． 四、实验内容 实验内容主要包括实验的实施方案、步骤、实验数据准备、实验的算法以及可能用到的仪器设备． 五、程序流程图 画出程序实现过程的流程图，以便更好的对程序执行的过程有清楚的认识，在程序调试过程中更容易发现问题． 六、实验结果 实验结果应包括实验的原始数据、中间结果及实验的最终结果，复杂的结果可以用表格

形式列出，较为简单的结果可以与实验结果分析合并出现． 七、实验结果分析 实验结果分析包括对对算法的理解与分析、改进与建议． 数值实验报告范例 为了更好地做好数值实验并写出规范的数值实验报告，下面给出一简单范例供读者参考． 数值实验报告 小组名称: 小组成员(班号): 按贡献排序情况: 指导教师评语: 小组所得分数: 一、实验名称 误差传播与算法稳定性． 二、实验目的 1．理解数值计算稳定性的概念． 2．了解数值计算方法的必要性． 3．体会数值计算的收敛性与收敛速度． 三、实验内容 计算dx x x I n n ? += 1 10 ，1,2,,10n = ． 四、算法描述 由 dx x x I n n ? += 1 10 ，知 dx x x I n n ?+=--101110，则

《数值计算方法》上机实验报告

《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程，若已知根的一个近似值，将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项，有 ，fxfxfxxx 0 ()()()，，, kkk 右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的 **根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0,

解出 fX 0 *k XX,, k' fx 0 k 水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:，见, ,3,输入变量、输出变量说明: X输入变量:迭代初值，迭代精度，迭代最大次数，\0 输出变量:当前迭代次数，当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx，，01* fx ()0 XX,,，?10 kk, ,1，kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志

,3,输入变量、输出变量说明: 结束 例：导出计算的牛顿迭代公式，并il ?算。（课本P39例2-16） 115cc （0）, 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组，1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝 aakkkkkk 对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行（包括常ekk 数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结 ,2,计算机程序框图:，见下页, 输入变量:系数矩阵元素，常向量元素baiji 输出变量:解向量元素bbb,,12n

(完整版)哈工大-数值分析上机实验报告

实验报告一 题目：非线性方程求解 摘要：非线性方程的解析解通常很难给出，因此线性方程的数值解法就尤为重要。本实验采用两种常见的求解方法二分法和Newton法及改进的Newton法。 前言：（目的和意义） 掌握二分法与Newton法的基本原理和应用。 数学原理： 对于一个非线性方程的数值解法很多。在此介绍两种最常见的方法：二分法和Newton法。 对于二分法，其数学实质就是说对于给定的待求解的方程f(x)，其在[a,b]上连续，f(a)f(b)<0，且f(x)在[a,b]内仅有一个实根x*，取区间中点c，若，则c恰为其根，否则根据f(a)f(c)<0是否成立判断根在区间[a,c]和[c,b]中的哪一个，从而得出新区间，仍称为[a,b]。重复运行计算，直至满足精度为止。这就是二分法的计算思想。

Newton法通常预先要给出一个猜测初值x0，然后根据其迭代公式 产生逼近解x*的迭代数列{x k}，这就是Newton法的思想。当x0接近x*时收敛很快，但是当x0选择不好时，可能会发散，因此初值的选取很重要。另外，若将该迭代公式改进为 其中r为要求的方程的根的重数，这就是改进的Newton法，当求解已知重数的方程的根时，在同种条件下其收敛速度要比Newton法快的多。 程序设计： 本实验采用Matlab的M文件编写。其中待求解的方程写成function的方式，如下 function y=f(x); y=-x*x-sin(x); 写成如上形式即可，下面给出主程序。 二分法源程序： clear %%%给定求解区间 b=1.5; a=0;

%%%误差 R=1; k=0;%迭代次数初值 while (R>5e-6) ; c=(a+b)/2; if f12(a)*f12(c)>0; a=c; else b=c; end R=b-a;%求出误差 k=k+1; end x=c%给出解 Newton法及改进的Newton法源程序：clear %%%% 输入函数 f=input('请输入需要求解函数>>','s') %%%求解f(x)的导数 df=diff(f);

计算方法实验报告 拟合

南京信息工程大学实验（实习）报告 一、实验目的: 用最小二乘法将给定的十个点拟合成三次多项式。 二、实验步骤： 用matlab编制以函数为基的多项式最小二乘拟合程序，并用于对下列数据作三次多项式最小二乘拟合(取权函数wi=1) x -2.5 -2 -1.5 -1 -0.5 0 0.5 1 1.5 2 y -2.30 -1 -0.14 -0.25 0.61 1.03 1.75 2.75 4.42 6.94 给定直线方程为：y=1/4*x3+1/2*x2+x+1 三、实验结论： 最小二乘法：通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据，并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。 一般地。当测量数据的散布图无明显的规律时，习惯上取n次代数多项式。 程序运行结果为： a = 0.9731 1.1023 0.4862 0.2238 即拟合的三次方程为：y=0.9731+1.1023x+0.4862*x2+0.2238*x3

-2.5 -2-1.5-1-0.5 00.51 1.52 2.5 -4-20246 81012 x 轴 y 轴 拟合图 离散点 y=a(1)+a(2)*x+a(3)*x.2+a(4)*x.3 结论： 一般情况下，拟合函数使得所有的残差为零是不可能的。由图形可以看出最小二乘解决了残差的正负相互抵消的问题，使得拟合函数更加密合实验数据。 优点：曲线拟合是使拟合函数和一系列的离散点与观测值的偏差平方和达到最小。 缺点：由于计算方法简单，若要保证数据的精确度，需要大量的数据代入计算。

西安交通大学计算方法B大作业

计算方法上机报告 姓名： 学号： 班级：

目录 题目一------------------------------------------------------------------------------------------ - 4 - 1.1题目内容 ---------------------------------------------------------------------------- - 4 - 1.2算法思想 ---------------------------------------------------------------------------- - 4 - 1.3Matlab源程序----------------------------------------------------------------------- - 5 - 1.4计算结果及总结 ------------------------------------------------------------------- - 5 - 题目二------------------------------------------------------------------------------------------ - 7 - 2.1题目内容 ---------------------------------------------------------------------------- - 7 - 2.2算法思想 ---------------------------------------------------------------------------- - 7 - 2.3 Matlab源程序---------------------------------------------------------------------- - 8 - 2.4计算结果及总结 ------------------------------------------------------------------- - 9 - 题目三----------------------------------------------------------------------------------------- - 11 - 3.1题目内容 --------------------------------------------------------------------------- - 11 - 3.2算法思想 --------------------------------------------------------------------------- - 11 - 3.3Matlab源程序---------------------------------------------------------------------- - 13 - 3.4计算结果及总结 ------------------------------------------------------------------ - 14 - 题目四----------------------------------------------------------------------------------------- - 15 - 4.1题目内容 --------------------------------------------------------------------------- - 15 - 4.2算法思想 --------------------------------------------------------------------------- - 15 - 4.3Matlab源程序---------------------------------------------------------------------- - 15 - 4.4计算结果及总结 ------------------------------------------------------------------ - 16 - 题目五----------------------------------------------------------------------------------------- - 18 -

数值计算方法上机实习题

数值计算方法上机实习题 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+ -=-，从I 0=0.1824, 0=0.1823I 出发，计算20I ； （2） 20=0I ，20=10000I , 用n I I n n 51 5111+- =--，计算0I ； （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 答：第一个算法可得出 e 0=|I 0?I 0 ?| e n =|I n ?I n ?|=5n |e 0| 易知第一个算法每一步计算都把误差放大了5倍，n 次计算后更是放大了5n 倍，可靠性低。 第二个算法可得出 e n =|I n ?I n ?| e 0=(15 )n |e n | 可以看出第二个算法每一步计算就把误差缩小5倍，n 次后缩小了5n 倍，可靠性高。

2． 求方程0210=-+x e x 的近似根，要求41105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 计算根与步数程序： fplot(@(x) exp(x)+10*x-2,[0,1]); grid on; syms x; f=exp(x)+10*x-2; [root,n]=EFF3(f,0,1); fprintf('root=%6.8f ,n=%d \n',root,n); 计算结果显示： root=0.09057617 ,n=11 （2） 取初值00=x ，并用迭代10 21 x k e x -=+；

（3） 加速迭代的结果; （4） 取初值00 x ，并用牛顿迭代法；

西交大计算方法上机报告

计算方法（B）实验报告 姓名： 学号： 学院： 专业：

实验一 三对角方程组Tx f =的求解 一、 实验目的 掌握三对角方程组Tx f =求解的方法。 二、 实验内容 求三对角方程组Tx f =的解，其中： 4 -1 -1 4 -1 -1 4 1 -1 4T ????????=?? ?? ???? ， 3223f ?? ? ? ?= ? ? ??? 三、 算法组织 设系数矩阵为三对角矩阵 11222333111 b c a b c a b c a b c b n n n n T ---???????? =?????? ?????? 则方程组Tx f =称为三对角方程组。 设矩阵T 非奇异，T 可分解为T=LU ，其中L 为下三角矩阵，U 为单位上三角矩阵，记 1 1 212 313 1 1 1111 ,11n n n n n r l r l r L U l r l μμμμμ---???? ? ? ? ? ? ?== ? ? ? ? ? ? ? ? ? ?? ? ? ? 可先依次求出,L U 中的元素后，令Ux y =，先求解下三角方程组Ly f =得出 y ，再求解上三角方程组Ux y =。 追赶法的算法组织如下： 1.输入三对角矩阵T 和右端向量f ；

2.将Tx f =压缩为四个一维数组{}{}{}{}i i i i a b c d 、、、，{}{}{}i i i a b c 、、是T 的三对角线性方程组的三个对角，{}i d 是右端向量。将分解矩阵压缩为三个一维数组 {}{}{}i i i l r μ、、。 3.对T 做Crout 分解（也可以用Doolittle 分解）导出追赶法的计算步骤如下： 1111,b r c μ== for 2i n = 111, , ,i i i i i i i i i i i i i l a b a r r c y d l y μμ---==-==- end 4.回代求解x /n n n x y μ= for 11i n =- 1()/i i i i i x y c x μ+=- end 5. 停止，输出结果。 四、 MATLAB 程序 MATLAB 程序见附件1. 五、 结果及分析 实验结果为： (1.0000 1.0000 1.0000 1.0000)T x =

计算方法B上机报告

计算方法B 上机报告 第1题 某通信公司在一次施工中，需要在水面宽度为20米的河沟底部沿直线走向铺设一条沟底光缆。在铺设光缆之前需要对沟底的地形进行初步探测，从而估计所需光缆的长度，为工程预算提供依据。已探测到一组等分点位置的深度数据(单位:米)如下表所示： (1)请用合适的曲线拟合所测数据点； (2)估算所需光缆长度的近似值，并作出铺设河底光缆的曲线图； 问题分析和算法思想： 本题的主要目的是对21个测量数据进行拟合，同时对拟合曲线进行线积分即可得到河底光缆长度的近似值，可以用的插值方法很多：多项式插值、Lagrange 插值、Newton 插值、三次样条插值等。由于数值点较多时，采用高次多项式插值将产生很大的误差，用拉格朗日插值多项式会出现龙格现象。故为了将所有的数据点都用上，且题中光缆为柔性，可光滑铺设于水底，鉴于此特性，采用三次样条插值的方法较为合适。 计算光缆长度近似值，只需将每两点之间的距离算出，然后依次相加，所得的折线长度，即为光缆长度的近似值。 光缆长度计算公式： 19 1 k k k l +===∑? ? ? 算法结构： 三次样条算法结构见《计算方法教程》P110。 源程序： clear;clc; x=0:20;

y=[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]; d=y; plot(x,y,'k.','markersize',15) hold on %%%计算二阶差商 for k=1:2 for i=21:-1:(k+1) d(i)=(d(i)-d(i-1))/(x(i)-x(i-k)); end end %%%假定d的边界条件，采用自然三次样条 for i=2:20 d(i)=6*d(i+1); end d(1)=0; d(21)=0; %%%追赶法求解带状矩阵的m值 a=0.5*ones(1,21); b=2*ones(1,21); c=0.5*ones(1,21); a(1)=0;c(21)=0; u=ones(1,21); u(1)=b(1); r=c; yy(1)=d(1); %%%追的过程 for k=2:21 l(k)=a(k)/u(k-1); u(k)=b(k)-l(k)*r(k-1); yy(k)=d(k)-l(k)*yy(k-1); end %%%赶的过程 m(21)=yy(21)/u(21); for k=20:-1:1 m(k)=(yy(k)-r(k)*m(k+1))/u(k); end %%%利用插值点画出拟合曲线 k=1; nn=100; xx=linspace(0,20,nn); l=0; for j=1:nn for i=2:20 if xx(j)<=x(i) k=i;