计算方法A 上机大作业
张晓璐 硕4011班 学号：3114009097
1. 共轭梯度法求解线性方程组
算法原理：由定理3.4.1可知系数矩阵A 是对称正定矩阵的线性方程组Ax=b
的解与求解二次函数1()2T T
f x x Ax b x =- 极小点具有等价性，所以可以利用
共轭梯度法求解1()2
T T
f x x Ax b x =-的极小点来达到求解Ax=b 的目的。
共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式：
(1)()()k k k k x x d α+=+
产生的迭代序列(1)(2)(3)x x x ，，，... 在无舍入误差假定下，最多经过n 次迭代，就可求得()f x 的最小值，也就是方程Ax=b 的解。 首先导出最佳步长k α的计算式。
假设迭代点()k x 和搜索方向()k d 已经给定，便可以通过()()()()k k f x d φαα=+ 的极小化
()()min ()()k k f x d φαα=+
来求得，根据多元复合函数的求导法则得:
()()()'()()k k T k f x d d φαα=?+
令'()0φα=，得到:
()()
()()k T k k k T k r d d Ad
α= ，其中()()k k r b Ax =-
然后确定搜索方向()k d 。给定初始向量(0)x 后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0)(0)(0)(0)()d r f x b Ax ==-?=- 。令
(1)(0)00x x d α=+
其中(0)(0)
0(0)(0)T T r d d Ad
α=。第二次迭代时，从(1)x 出发的搜索方向不再取(1)r ，而
是选取(1)
(1)(0)0d
r d β=+，使得(1)d 与(0)d 是关于矩阵A 的共轭向量，由此可求
得参数0β:
(1)(0)
0(0)(0)T T r Ad d Ad
β=-
然后从(1)x 出发，沿(1)d 进行搜索得到(2)
(1)(1)1x x d α=+
设已经求出(1)
()()k k k k x x d α+=+，计算(1)(1)k k r b Ax ++=-。
令(1)
(1)()k k k k d
r d β++=+，选取k β，使得(1)k d +和()k d 是关于A 的共轭向量，可
得:
(1)()
()()k T k k k T k r Ad d Ad
β+=-
具体编程计算过程如下：
（1） 给定初始近似向量(0)x 以及精度0ε> ； （2） 计算(0)(0)r b Ax =-，取(0)(0)d r =； （3） For k=0 to n-1 do
（i ）()()
0()()k T k k T k r d d Ad
α=；
（ii ）(1)
()()k k k k x
x d α+=+；
（iii ）(1)(1)k k r b Ax ++=-；
（iv ）若(1)
k r ε+≤或k=n-1，则输出近似解(1)k x + ，停止；否则，转（v ）；
（v ）2
(1)22()
2
k k k r r
β+=
；
（vi ）(1)
(1)()k k k k d
r d β++=+；
End do
程序框图：
程序使用说明：本共轭梯度法求解线性方程的程序是一个函数Gongetidu2(A,b,x0,epsa)，在求解线性方程组Ax=b（A是对称正定矩阵）的时候，首先给定初始向量x0和误差epsa，然后直接调用该函数Gongetidu2(A,b,x0,epsa)即可，其中A,b,x0和epsa都是可变的，虽然该函数是通用的，但是对于矩阵A和向量b的输入，需要使用者根据A和b的特点自行输入。
算例3.2计算结果：
对题113页3.2，首先编程将矩阵A,b,x0,epsa读入系统，然后再调用函数Gongetidu2(A,b,x0,epsa)；
程序如下：
clear A b x0 %程序运行前首先清除A，b和X0的数值，以免影响计算
clc
n=input('请输入对称正定矩阵A的阶数n=');
epsa=input('请输入误差=');
for k=1:(n-1) %读取题目3.2的矩阵A
A(k,k)=-2;
A(k,k+1)=1;
A(k+1,k)=1;
end
A(n,n)=-2; A;
b(1,1)=-1; %读取题目3.2的矩阵b b(n,1)=-1; b;
x0(1,1)=input('假定初始向量各元素相同，均为：'); %给题目3.2的迭代过程赋初值
for kk=2:n
x0(kk,1)=x0(1,1); end
x=Gongetidu2(A,b,x0,epsa) %调用共轭梯度法求线性方程组的函数
function x=Gongetidu2(A,b,x0,epsa) n=size(A,1); x=x0; r=b-A*x; d=r;
for k=0:(n-1)
alpha=(r'*r)/(d'*A*d); x=x+alpha*d; r2=b-A*x;
if ((norm(r2)<=epsa)|(k==n-1)) x; break; end
beta=norm(r2)^2/norm(r)^2; d=r2+beta*d; r=r2; end
计算结果：当n=100时，x=[1;1;1;1;1;1;1;…..1;1;1;1]
当n=200时，x=[1;1;1;1;1;1;1;…..1;1;1;1;1;1;1]
当n=400时，x=[1;1;1;1;1;1;1;…..1;1;1;1;1;1;1;1;1;1] 2. 三次样条差值
算法原理（三次样条插值函数的导出）： (i).导出在子区间[]1,i i x x - 上的S(x)的表达式
由于S(x)的二阶导数连续，设S(x)再节点i x 处的二阶导数值S ’’(xi)=Mi ，其中Mi 为未知的待定参数。由S(x)是分段的三次多项式知，S ’’(x)是分段线性函数，S ’’(x)在子区间[]1,i i x x -上可表示为
1
111
111''(),i i i i
i i i i i i i i i i
i i
x x x x S x M M x x x x x x x x M M x x x h h ---------=
+----=+≤≤
其中hi=xi-x(i-1),对上式两次积分得到
()()3
3
11
11()66()(),i
i i i i
i
i i i i i i
x x x x S x M M h h b x x c x x x x x ------=+
+
-+-≤≤
由插值条件11(),()i i i i S x y S x y --== 得到
22
11()/,()/66
i i i i i i i i i i h h b y M h c y M h --=-=-
将i i b c 和 代入()S x 可得
()()3
3
2
11
112
11()()/()
666
()/(),6
i
i i i i i i i i i i
i
i i i i i i x x x x h S x M M y M h x x h h h
y M h x x x x x --------=+
+--+-
-≤≤
(ii).建立参数i M 的方程组 对
S(x)求导可得
()()2
2
11
11
1'()()/22,6
i
i i i i i i
i
i
i i i i i
x x x x S x M M y y h h h M M h x x x -------=-+
+---
≤≤
上式中令i x x = 得S(x)在xi 处的左导数'
()i S x - ，令1i x x -=得到右导数'
()i S x +，因为S(x)在内节点xi 处一阶导数连续，所以'
'
()()i i S x S x -+=，进一步推导可得
112,1,2,3,...,1i i i i i i M M M d i n μλ-+++==-
其中
1
i
i i i h h h μ+=
+，
1
1
1i i i i i h h h λμ++=
=-+，
1111116
()6[,,],1,2,...,1i i i i i i i i i i i i
y y y y d f x x x i n h h h h +--+++--=
-==-+
上式为三弯矩方程组，因为三弯矩方程组只有n-1个方程，不能确定n+1个未知
量Mi ，所以需要再增加两个方程，由边界条件确定。
第一种边界条件:此时已知''()''()f a f b 和 .不妨取0''()M f a =，''()n M f b =，这时三弯矩方程组化为：
1121101112
111222,3,...,22i i i i i i
n i n n n n M M d M M M M d i n M M d M
λμμλμλ-+-----+=-??
++==-??+=-? 以上方程组系数矩阵式严格三对角占优矩阵，可用追赶法求解。求出
(1,2,...,1i M i n =- 后，代入S(x)可得三次样条插值函数的数学表达式。 第二种边界条件：已知'()'()f a f b 和。记0''()''()n y f a y f b ==,，则有
00'()''()'n n S x y S x y +-==, 所以：
1011101011',',3663n n n n n n n n
y y h h y y h h
M M y M M y h h -----
-+=++= 即
0102M M d +=
12n n n M M d -+=
其中
10
0011
16(')6
(')
n n n n n n
y y d y h h y y d y h h --=
--=-
所以得到第二种边界条件下的三弯矩方程组：
010111
2,
2,1,2,3,...,1,2i i i i i i n n n M M d M M M d i n M M d
μλ-+-+=??
++==-??+=?
该方程组系数矩阵是严格三对角占优矩阵，可用追赶法求解，具体追赶法的求解
过程见《数值分析》教材。
第三种边界条件：周期型边界条件.已知()y f x = 是以0n T b a x x =-=- 为周期的周期函数，则由周期性可知，01101111,,,,n n n n n y y y y M M M M h h +++===== ，
这时，可以将点n x 看成内点，则方程组对i=n 也成立，既有
112n n n n n n M M M d μλ-+++= ，
也即
112n n n n n M M M d λμ-++=，
其中
11116
()n n n n n n
y y y y d h h h h ---=
-+ 于是三弯矩方程组化为
11211111
12,
2,2,3,4,....,1,2.
n i i i i i i n n n n n M M M d M M M d i n M M M d λμμλλμ-+-++=??
++==-??++=? 该方程组可用LU 分解法求解，具体追赶法的求解过程见《数值分析》教材。
程序框图如下：
程序使用说明：因为该程序需要计算三种不同边界条件的三弯矩方程组，所以首先定义追赶法和LU 分解法求解线性方程组的函数，ZhuiGanFa(A,b)和LU_fenjieqiuxianxingfangcheng(A,b)，然后在确定了边界条件类型之后，构造关于M 的三弯矩方程AM=b 的A 和b 矩阵，然后分别代入函数ZhuiGanFa(A,b)和LU_fenjieqiuxianxingfangcheng(A,b)便可求得M ，然后根据
332111121
1()()()()666(),6i i i i
i i i i i i i
i
i i i i i i
x x x x h x x S x M M y M h h h h x x y M x x x h ---------=++--+-
≤≤
便可得到，在1i i x x x -≤≤ 上的三次样条插值函数()S x ，进而得到整个区间上的三次样条差值函数()S x 。
算例计算结果： 141页的计算实习
当为第一类边界条件时：
当-1= 当为第二类边界条件时： 当-1= x y 三次样条差值函数S10(x)和原函数f(x)的对比 当0= 3. 龙贝格积分法：对于复化梯形求积公式，取积分近似值 2221 ()()41 n n n n I f T T T T ≈+-=- 对复化辛普森求积公式 4(4) ()(),2880 n b a I f S h f a b ηη--≈-≤≤ 4(4) 211()()(),28802 n b a h I f S f a b ηη--≈-≤≤ 因为(4)(4)1()()f f ηη=，所以上述两式相除得 2()16()n n I f S I f S -≈- 所以， 2222 1 ()()41 n n n n I f S S S S ≈+ -=- 同理， 22231 ()()41 n n n n I f C C C C ≈+ -=- 对2n T ，2n S 和2n C 分析可得 22222 2231()41 1()411()41n n n n n n n n n n n n S T T T C S S S R C C C ? =+-?-? ? =+-?-? ? =+-?-? x y 三次样条差值函数S10(x)和原函数f(x)的对比 龙贝格积分算法如下： 11111111211221 2 22222222322 22[()()], 2 1((21)),0,1,2...,2221(),41 1(), 0,1,2..., 411(),41k k k k k k k k k k k k k k k k k i b a T f a f b b a b a T T f a i k S T T T C S S S k R C C C +++++++++=-=+--=++-=?=+-?-? ? =+-=?-? ? =+-?-? ∑ 如果122,k k R R ε+-< 则取12[]k I f R +≈，否则，继续计算直到满足条件。 程序框图： 程序使用说明:运行本程序的时候，直接按照提示输入所求积分的原函数()f x （比如11x + ），然后按提示依次输入积分下限a ，积分上限b 和积分 精度1eps ，然后程序便可计算出原函数()f x 在[,]a b 之间的积分数值。 算例计算结果： 6.2 101 0.69266061dx x =+? 120ln(1) 0.27291381x dx x +=+? 10ln(1) 0.8222677x dx x +=? /20sin() 1.3714136x dx x π=? 4. 四阶龙格-库塔法求解常微分方程的初值问题 算法原理：用标准4级4阶R-K 法求解一阶常微分方程的算法公式为： 1 12341213 2431(22)6 (,)11(,)2211(,)22 (,) i i i i i i i i i i y y K K K K K hf x y K hf x h y K K hf x h y K K hf x h y K +? =++++?? =??? =++?? ? =++??=++?? 求解m 阶常微分方程 ()(1)(1)(1) 000(,,','',...,);,(),'()',...,()m m m m y f x y y y y a x b y a y y a y y a y ---?=≤≤??===?? 将其转化为一阶微分方程组来求解。为此，引进新的变量1y y = , 2'y y =,…, (1)m m y y -=,即可将m 阶微分方程，转化为如下的一阶微分方程组： 12 2 31 12,...,(1)10200 ''... ''(,,)(),()',...,()m m m m m m y y y y y y y f x y y y y a y y a y y a y --=??=????=??=?===?? 程序框图： 程序使用说明：运行该程序，按照提示依次输入常微分方程的阶数m ，x 的下 限a 和上限b ，步长h 以及y 的最高阶导数表达式()(1)(2)...()m m m y ay by f x --=+++ ， 特别需要注意的是：由于程序设定的问题，一定要将公式中的()k y 用(1,1)y k +代替，例如所求解的初值问题为 2''2'2sin 01x y y y e x x -+=≤≤， ，那么本程序所需要输入的公式为：2*(2,1)2*(1,1)exp(2*)*sin()y y x x -+ 。且()f x 代表含x 的表达式。接着依次输入已知条件()y a ，'()y a ，''()y a ，…，即可，最终可以得出该问题的 解。 算例计算结果： 9.2题结果 (0)-1.000000000000000(0.05)-0.847436458333333(0.1)-0.689482936121419(0.15)-0.525724908067489(0.2)-0.355719511320969(0.25)-0.178993729355764(0.3)0.004957535888930(0.35)0.196673510288970y y y y y y y y ========(0.4)0.396729719312157(0.45)0.605740292781225(0.5)0.824360410550626(0.55) 1.053288897488131(0.6) 1.293270976642268(0.65) 1.545101189994835(0.7) 1.809626496745704(0.75) 2.087749559656611(0.y y y y y y y y y ========8) 2.380432230593372(0.85) 2.688699247053874(0.9) 3.013642152154253(0.95) 3.356423451269989(1) 3.718281019294475 y y y y ===== x y 近似解和解析解y=x*e x +2*x-1的比较 西安交大统计学考试试卷 一、单项选择题（每小题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、“性别”是品质标志。（对） 第一章 误差 1 问3.142，3.141，7 22分别作为π的近似值各具有几位有效数字？ 分析 利用有效数字的概念可直接得出。 解 π=3.141 592 65… 记x 1=3.142，x 2=3.141，x 3=7 22. 由π- x 1=3.141 59…-3.142=-0.000 40…知 34111 10||1022 x π--?<-≤? 因而x 1具有4位有效数字。 由π- x 2=3.141 59…-3.141=-0.000 59…知 223102 1||1021--?≤- 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; 计算方法上机报告 姓名： 学号： 班级： 上课班级： 说明： 本次上机实验使用的编程语言是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.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量 （1）二分法 （局部，大图不太看得清，故后面两小题都用局部截图） （2）迭代法 （3）牛顿法 顺序消元法 #include for(k=p+1;k 《数值计算方法》上机实验报告华北电力大学 实验名称数值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 计算方法第四次上机报告 2.用欧拉方法解初值 y’=10x(1-y) 0<=x<=1 Y(0)=0 取步长h=0.1,保留5位有效数字，并与准确解相比较 分析：该题目考察欧拉方法解初值问题 程序如下： function Heun(a,b,y0,n) h=(b-a)/n;x=a:h:b; y=y0*ones(1,n+1); for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1)); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc); end for k=1:n+1 fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end function z=f(xx,yy) z=10*xx*(1-yy); 运行结果： >> Heun(0,1,0,10) x[0]=0.000000 y[0]=0.000000 x[1]=0.100000 y[1]=0.050000 x[2]=0.200000 y[2]=0.183000 x[3]=0.300000 y[3]=0.362740 x[4]=0.400000 y[4]=0.547545 x[5]=0.500000 y[5]=0.705905 x[6]=0.600000 y[6]=0.823543 x[7]=0.700000 y[7]=0.901184 x[8]=0.800000 y[8]=0.947627 x[9]=0.900000 y[9]=0.973290 x[10]=1.000000 y[10]=0.986645 >> 分析： 该结果与准确结果相比比较接近，但是有一定的误差。 6.用四阶龙格—库塔公式解第三题中的初值问题，取步长h=0.2，保留五位有效数字。 题目目的分析： 该题考查四阶龙格-库塔方法和改进欧拉方法求解精确度问题。 程序： 改进欧拉法： function Heun(a,b,y0,n) h=(b-a)/n;x=a:h:b; y=y0*ones(1,n+1); for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1)); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc); end for k=1:n+1 fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end 计算方法 一、填空题 1.假定x ≤1，用泰勒多项式?+??+++=! !212n x x x e n x ,计算e x 的值，若要求截断误差不超过0.005，则n=_5___ 2. 解 方 程 03432 3=-+x － x x 的牛顿迭代公式 )463/()343(121121311+--+--=------k k k k k k k x x x x x x x 3.一阶常微分方程初值问题 ?????= ='y x y y x f y 0 0)() ,(,其改进的欧拉方法格式为)],(),([21 1 1 y x y x y y i i i i i i f f h +++++= 4.解三对角线方程组的计算方法称为追赶法或回代法 5. 数值求解初值问题的四阶龙格——库塔公式的局部截断误差为o(h 5 ) 6.在ALGOL 中，简单算术表达式y x 3 + 的写法为x+y ↑3 7.循环语句分为离散型循环,步长型循环,当型循环. 8.函数)(x f 在[a,b]上的一次（线性）插值函数= )(x l )()(b f a b a x a f b a b x --+-- 9.在实际进行插值时插值时，将插值范围分为若干段，然后在每个分段上使用低阶插值————如线性插值和抛物插值，这就是所谓分段插值法 10、数值计算中，误差主要来源于模型误差、观测误差、截断误差和舍入误差。 11、电子计算机的结构大体上可分为输入设备 、 存储器、运算器、控制器、 输出设备 五个主要部分。 12、算式2 cos sin 2x x x +在ALGOL 中写为))2cos()(sin(2↑+↑x x x 。 13、ALGOL 算法语言的基本符号分为 字母 、 数字 、 逻辑值、 定义符四大 计算方法实验报告 班级： 学号： 姓名： 成绩： 1 舍入误差及稳定性 一、实验目的 （1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； （2）通过上机计算，了解舍入误差所引起的数值不稳定性 二、实验内容 1、用两种不同的顺序计算10000 21n n -=∑，分析其误差的变化 2、已知连分数() 1 01223//(.../)n n a f b b a b a a b =+ +++，利用下面的算法计算f ： 1 1 ,i n n i i i a d b d b d ++==+ (1,2,...,0 i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分 1 041 n n x y dx x =+? (0,1,...,1 n = 4、设2 2 11N N j S j == -∑ ，已知其精确值为1311221N N ?? -- ?+?? （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序 （3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数 三、实验步骤、程序设计、实验结果及分析 1、用两种不同的顺序计算10000 2 1n n -=∑，分析其误差的变化 （1）实验步骤： 分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计： a.顺序计算 #include 习题一 1. 什么叫数值方法？数值方法的基本思想及其优劣的评价标准如何？ 数值方法是利用计算机求解数学问题近似解的方法 x max x i , x ( x 1 , x 2 , x n ) T R n 及 A n R n n . 2. 试证明 max a ij , A ( a ij ) 1 i n 1 i n 1 j 证明： （ 1）令 x r max x i 1 i n n p 1/ p n x i p 1/ p n x r p 1/ p 1/ p x lim( x i lim x r [ ( ] lim x r [ lim x r ) ) ( ) ] x r n p i 1 p i 1 x r p i 1 x r p 即 x x r n p 1/ p n p 1/ p 又 lim( lim( x r x i ) x r ) p i 1 p i 1 即 x x r x x r ⑵ 设 x (x 1,... x n ) 0 ，不妨设 A 0 ， n n n n 令 max a ij Ax max a ij x j max a ij x j max x i max a ij x 1 i n j 1 1 i n j 1 1 i n j 1 1 i n 1 i n j 1 即对任意非零 x R n ，有 Ax x 下面证明存在向量 x 0 0 ，使得 Ax 0 ， x 0 n ( x 1,... x n )T 。其中 x j 设 j a i 0 j ，取向量 x 0 sign(a i 0 j )( j 1,2,..., n) 。 1 n n 显然 x 0 1 且 Ax 0 任意分量为 a i 0 j x j a i 0 j ， i 1 i 1 n n 故有 Ax 0 max a ij x j a i 0 j 即证。 i i 1 j 1 3. 古代数学家祖冲之曾以 355 作为圆周率的近似值，问此近似值具有多少位有效数字？ 113 解： x 325 &0.314159292 101 133 x x 355 0.266 10 6 0.5 101 7 该近似值具有 7 为有效数字。 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 计算方法（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 = 第一次&第二次上机作业 上机作业： 1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5 给出执行结果，并简要分析一下产生现象的原因。 解：执行结果如下： 在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。 2.（课本181页第一题） 解：（1）n=0时，积分得I0=ln6-ln5，编写如下图代码 从以上代码显示的结果可以看出，I 20的近似值为0.7465 （2）I I =∫I I 5+I 10dx,可得∫I I 610dx ≤∫I I 5+I 10dx ≤∫I I 510dx,得 16(I +1)≤I I ≤15(I +1),则有1126≤I 20≤1105, 取I 20=1 105 ,以此逆序估算I 0。代码段及结果如下图所示 （3）从I20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设I I表示I I的近似值，I I-I I=(?5)I(I0?I0)(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。 2.（课本181页第二题） （1）上机代码如图所示 求得近似根为0.09058 （2）上机代码如图所示 得近似根为0.09064； （3）牛顿法上机代码如下 计算所得近似解为0.09091 第三次上机作业上机作业181页第四题 线性方程组为 [1.13483.8326 0.53011.7875 1.16513.4017 2.53301.5435 3.4129 4.9317 1.23714.9998 8.76431.3142 10.67210.0147 ][ I1 I2 I3 I4 ]=[ 9.5342 6.3941 18.4231 16.9237 ] （1）顺序消元法 A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435; 3.4129, 4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147]; b=[9.5342;6.3941;18.4231;16.9237]; 上机代码（函数部分）如下 function [b] = gaus( A,b )%用b返回方程组的解 B=[A,b]; n=length(b); RA=rank(A); RB=rank(B); 计算方法A 上机大作业 1. 共轭梯度法求解线性方程组 算法原理：由定理3.4.1可知系数矩阵A 是对称正定矩阵的线性方程组Ax=b 的解与求解二次函数1()2 T T f x x Ax b x =-极小点具有等价性，所以可以利用共轭梯度法求解1()2 T T f x x Ax b x = -的极小点来达到求解Ax=b 的目的。 共轭梯度法在形式上具有迭代法的特征，在给定初始值情况下，根据迭代公式： (1)()()k k k k x x d α+=+ 产生的迭代序列(1)(2)(3)x x x ，，，... 在无舍入误差假定下，最多经过n 次迭代，就可求得()f x 的最小值，也就是方程Ax=b 的解。 首先导出最佳步长k α的计算式。 假设迭代点()k x 和搜索方向()k d 已经给定，便可以通过()()()() k k f x d φαα=+的极小化 ()()min ()()k k f x d φαα=+ 来求得，根据多元复合函数的求导法则得: ()()()'()()k k T k f x d d φαα=?+ 令'()0φα=，得到: ()() ()()k T k k k T k r d d Ad α=，其中()()k k r b Ax =- 然后确定搜索方向()k d 。给定初始向量(0)x 后，由于负梯度方向是函数下降最快的方向，故第一次迭代取搜索方向(0) (0)(0)(0)()d r f x b Ax ==-?=-。令 (1)(0)00x x d α=+ 其中(0)(0)0(0)(0) T T r d d Ad α=。第二次迭代时，从(1) x 出发的搜索方向不再取(1)r ，而是选取(1) (1)(0)0d r d β=+，使得(1)d 与(0)d 是关于矩阵A 的共轭向量，由此可 求得参数0β: : 第一次作业 1.下列各数都是经过四舍五入得到的近似数，指出他们有几位有效数字，并写出绝对误差限。 9800107480.566.385031.01021.1*65*5*4*3*2*1=?=====x x x x x x 解： 1* 11011021.01021.1?==x ，有5位有效数字，绝对误差限为4-5-1105.0105.0?=?； 1-* 2 1031.0031.0?==x ，有2位有效数字，绝对误差限为3-2-1-105.0105.0?=?； 3* 3103856.06.385?==x ；有4位有效数字，绝对误差限为-14-3105.0105.0?=?； 2* 41056480.0480.56?==x ；有5位有效数字，绝对误差限为3-5-2105.0105.0?=?； ; 65* 5 107.0107?=?=x ；有1位有效数字，绝对误差限为51-6105.0105.0?=?； 4* 6 109800.09800?==x ；有4位有效数字，绝对误差限为5.0105.04-4=?。 2.要使20的近似值的相对误差限小于%1.0，要取几位有效数字 解：由于110447213595.047213595.420??=?=，设要取n 位有效数字，则根据 定理，有()()%1.01081 1021111 《计算方法》上机指导书 实验1 MATLAB 基本命令 1．掌握MATLAB 的程序设计 实验内容：对以下问题，编写M 文件。 (1) 生成一个5×5矩阵，编程求其最大值及其所处的位置。 (2) 编程求∑=20 1!n n 。 (3) 一球从100米高度自由落下，每次落地后反跳回原高度的一半，再落下。求它在 第10次落地时，共经过多少米？第10次反弹有多高？ 2．掌握MATLAB 的绘图命令 实验内容：对于自变量x 的取值属于[0，3π]，在同一图形窗口画出如下图形。 （1）1sin()cos()y x x =?； （2）21 2sin()cos()3 y x x =-； 实验2 插值方法与数值积分 1. 研究人口数据的插值与预测 实验内容：下表给出了从1940年到1990年的美国人口，用插值方法推测1930年、1965年、2010年人口的近似值。 美国人口数据 1930年美国的人口大约是123,203千人，你认为你得到的1965年和2010年的人口数字精确度如何？ 2．最小二乘法拟合经验公式 实验内容：某类疾病发病率为y ‰和年龄段x (每五年为一段，例如0~5岁为第一段，6~10岁为第二段……)之间有形如bx ae y =的经验关系，观测得到的数据表如下 （1）用最小二乘法确定模型bx ae y =中的参数a 和b 。 （2）利用MATLAB 画出离散数据及拟合函数bx ae y =图形。 3.复化求积公式 实验内容：对于定积分? +=1 02 4dx x x I 。 （1）分别取利用复化梯形公式计算，并与真值比较。再画出计算误差与n 之间的曲线。 （2）取[0，1]上的9个点，分别用复化梯形公式和复化辛普森公式计算，并比较精度。 实验一 矩阵的分解 一、实验目的 掌握矩阵的分解原理和一般方法，学会利用矩阵分解直接求解线性方程组。 二、实验内容 求矩阵() 2020 =ij A α?的T LDL 分解与Cholesky 分解，其中 ,min(,),ij i i j i j i j α=?=? ≠? 。 三、问题分析 1． Cholesky 分解 Cholesky 分解是针对被分解矩阵为对称正定的情况给出的。 分解步骤如下： 11g =1111/y b g =，1111i i g g α= 2i n = ； DO 2j n = jj g = IF 0jj g < STOP ，JUMP TO (5) DO 1i j n =+ 1 1j ij ik kj k ij jj g g g g α-=??- ? ? ?=∑ ji ij g g = 1 1j i ik k k i jj b g y y g -=??- ? ? ?=∑ END DO END DO 2． T LDL 分解 T LDL 分解是针对Cholesky 分解中的开平方运算进行的改进。 分解步骤如下： 11i i r α=，1111/i i r r r =，11y b = 1i n = DO 2i n = DO j i n = 1 1i ij ij ik kj k r l r α-=??=- ??? ∑ /ji ij ii l r r = 1 1i i i ik k k y b l b -=??=- ??? ∑ END DO END DO 四、matlab 求解 分别写出T LDL 分解和Cholesky 分解的函数程序gaijinsqrt.m 和.cholesky m ，调用格 式如下： 1. [index,x,r]=gaijinsqrt(A,b) 参数说明： A 和b 分别是线性代数方程组Ax =b 的系数矩阵和右端向量；输出x 为解向量。 [index,x,g]=Cholesky(A,b) 参数说明： A 和b 分别是线性代数方程组Ax =b 的系数矩阵和右端向量；输出x 为解向量。 然后写出主程序2homework .m 如下： %生成矩阵A A=zeros(20,20); for i=1:20 for j=1:20 if i~=j if i>j A(i,j)=j; else A(i,j)=i; end 练习题与答案 练习题一 练习题二 练习题三 练习题四 练习题五 练习题六 练习题七 练习题八 练习题答案 练习题一 一、是非题 1.*x=–1 2.0326作为x的近似值一定具有6位有效数字，且其误差限 ≤ 4 10 2 1 - ? 。() 2.对两个不同数的近似数，误差越小，有效数位越多。( ) 3.一个近似数的有效数位愈多，其相对误差限愈小。( ) 4.用 2 1 2 x - 近似表示cos x产生舍入误差。( ) 5. 3.14和 3.142作为π的近似值有效数字位数相同。 ( ) 二、填空题 1. 为了使计算 ()()2334912111y x x x =+ -+ ---的乘除法次数尽量少，应将该 表达式改写为 ； 2. * x =–0.003457是x 舍入得到的近似值，它有 位有效数字，误差限 为 ，相对误差限为 ； 3. 误差的来源是 ； 4. 截断误差为 ； 5. 设计算法应遵循的原则是 。 三、选择题 1．* x =–0.026900作为x 的近似值，它的有效数字位数为( ) 。 (A) 7； (B) 3； (C) 不能确定 (D) 5. 2．舍入误差是( )产生的误差。 (A) 只取有限位数 (B) 模型准确值与用数值方法求得的准确值 (C) 观察与测量 (D) 数学模型准确值与实际值 3．用 1+x 近似表示e x 所产生的误差是( )误差。 (A). 模型 (B). 观测 (C). 截断 (D). 舍入 4．用s *=21 g t 2表示自由落体运动距离与时间的关系式 (g 为重力加速度)，s t 是在 时间t 内的实际距离，则s t - s *是（ ）误差。 (A). 舍入 (B). 观测 (C). 模型 (D). 截断 5．1.41300作为2的近似值，有( )位有效数字。 (A) 3； (B) 4； (C) 5； (D) 6。 四、计算题 计算方法上机报告 姓名: 学号: 班级: 目录 题目一-------------------------------------------------------------------- 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- 5.1题目内容 --------------------------------------------------------- 18- 5.2算法思想 ---------------------------------------------------------- 18 5.3 Matlab 源程序 ------------------------------------------------------ 18 5.3.1 非压缩带状对角方程组----------------------------------- 18 - 5.3.2压缩带状对角方程组--------------------------------------- 20- 5.4实验结果及分析 ----------------------------------------------- 22-统计西安交大期末考试试题(含答案)
计算方法_习题第一、二章答案..
西安交通大学计算方法B上机试题
计算方法上机作业
计算方法上机题答案
《数值计算方法》上机实验报告
计算方法上机作业
计算方法试题库讲解
计算方法上机实习题大作业(实验报告).
(完整word版)计算方法习题集及答案.doc
(完整版)数值计算方法上机实习题答案
西交大计算方法上机报告
计算方法上机作业集合
西交计算方法A上机大作业
计算方法作业参考答案(不断更新)
计算方法作业2
西安交通大学计算方法A实验报告
计算方法练习题与答案
西安交通大学计算方法B大作业资料
相关主题
文本预览