数值计算方法大作业
题目利用数值计算方法求取基尼系数
姓名与学号
指导教师
年级与专业
所在学院
一、问题综述：
基尼系数（Gini coefficient），是20世纪初意大利学者科拉多·吉尼根据劳伦茨曲线所定义的判断收入分配公平程度的指标。是比例数值，在0和1之间。基尼指数（Gini index）是指基尼系数乘100倍作百分比表示。在民众收入中，如基尼系数最大为“1”，最小等于“0”。前者表示居民之间的收入分配绝对不平均（即所有收入都集中在一个人手里，其余的国民没有收入），而后者则表示居民之间的收入分配绝对平均，即人与人之间收入绝对平等，但这两种情况只出现在理论上；因此，基尼系数的实际数值只能介于0～1之间，基尼系数越小收入分配越平均，基尼系数越大收入分配越不平均。
设右图中的
实际收入分配曲线
（红线）和收入分
配绝对平等线（绿
线）之间的面积为
A，和收入分配绝
对不平等线（蓝
线）之间的面积为
B，则表示收入与
人口之间的比例的基尼系数为
A
A+B
。
如果A为零，即基尼系数为0，表示收入分配完全平等（红线和绿线重叠）；如果B为零，则系数为1，收入分配绝对不平等（红线和蓝线重叠）。该系数可在0和1之间取任何值。实际上，一般国家的收入分配，既不是完全平等，也不是完全不平等，而是在两者之间，劳伦茨曲线为一条凸向横轴的曲线。收入分配越趋向平等，劳伦茨曲线的弧度越小（斜度越倾向45度），基尼系数也越小；反之，收入分配越趋向不平等，劳伦茨曲线的弧度越大，那么基尼系数也越大。
基尼系数的调节需要国家通过财政政策进行国民收入的二次分配，例如对民众的财政公共服务支出和税收等，从而让收入均等化，令基尼系数缩小。
基尼系数由于给出了反映居民之间贫富差异程度的数量界线，可以较客观、直观地反映和监测居民之间的贫富差距，预报、预警和防止居民之间出现贫富两极分化。因此得到世界各国的广泛认同和普遍采用。
联合国有关组织规定：
●若低于0.2表示收入平均；
●0.2-0.3表示相对平均；
●0.3-0.4表示相对合理；
●0.4-0.5表示收入差距大；
●0.6以上表示收入差距悬殊。
2013年1月18日，中国国家统计局一次性公布了自2003年以来十年的全国基尼系数。大陆统计局局长马建堂称，按照国际新的统计口径，大陆居民收入的基尼系数，2003年是0.479，2004年是0.473，2005年为0.485，2006年为0.487，2007年为0.484，2008年为0.491，2009年为0.490，2010年为
0.481，2011年为0.477，到2012年的数据是0.474，为2005年以来最低水平，而自2008年起，基尼系数也在逐年下降。而此前西南财大调查数据显示，中国的2012年的基尼系数为0.61，但无论是民间统计的数据还是官方统计的数据，结果都遭到学术界质疑，仍具有争议性。
本文将根据网络上国家统计局的数据，利用上面给出的公式来计算我国从2002年以来的城镇居民基尼系数，并将计算出的数据与现有数据进行比较。
全球基尼系数
二、问题分析：
想要通过拟合曲线再积分求取基尼系数，需要先求取劳伦茨曲线。劳伦茨曲线是1905年由经济学家马克斯·劳伦茨所提出的表示收入分配的曲线。
在经济学里,劳伦兹曲线是在过往财富分配数据上建立的累积分布函数所对应的曲线,它通过变量y%的值来反映各项分配的比例。它经常被用来描述收入的分配情况，即以x%代表一部分(收入相似)家庭占整个社会家庭的比例，以y%代表该部分家庭的收入占整个社会收入的比例。该曲线也可用来描述社会资本的分配情况。在这些应用当中，经济学家经常把它用来衡量社会(主要指社会收入)是否公平。所利用的概率密度分布函数为：
本文中使用国家统计局网站上面的数据，先采用拟合得方法来得到劳伦茨曲线，然后再用积分的方法来计算我国城镇居民的基尼系数。并分别使用多种方法进行计算，例如直接用求得的数据进行插值积分来计算基尼系数，比较不同，再和上文中所提到的国家统计局统计出的基尼系数进行比较，进行进一步的分析。
三、问题解决
首先在国家统计局的网站上面找到了按收入等级分城镇居民家庭基本情况——城镇居民平均每人全部年收入。列举数据表如下（困难户是包含在最低收入户中的群体）：
表1.2007-2012年城镇居民平均每人全部年收入
表2.2002-2006年城镇居民平均每人全部年收入
统计局的数据将所有城镇居民按收入划分为八个部分，按照劳伦茨曲线的定义，我们对每年的数据进行分析都可以得到七个数据点。通过这七个数据点再加上固有的（0,0），（1,1）点，每年共9个点，就是我们需要的数据。
数据处理方法：
从低到高算起，拿2012年的数据为例，对于城镇居民困难户（5%），其收入占总收入的比例为：
7520.90?5%
=1.39%
26959.00?100%
如果令第一个点为固定点（0,0），则由上面的结果可以求得第二个点为：（0.0139,0.05）。
同样，对于城镇居民最低收入户(10%)，其收入占总收入的比例为：
9209.50?10%
=3.42%
26959.00?100%
由于困难户是包含在最低收入户中的，所以不用进行累加，可直接得到第三个点：（0.0342,0.10）。
之后，对于城镇居民较低收入户(10%)，其收入占总收入的比例为：
13724.70?10%
=5.09%
26959.00?100%
由于与之前不再具有包含关系，所以需要与之前的计算结果进行累加，可得到第四个点：（0.0851,0.20）。之后的计算方式以此类推，此处不再赘述。
以此种方法对以上数据处理后，可得到如下结果：
表3.2007-2012年数据处理结果
表4.2002-2006年数据处理结果
1.最小二乘法拟合
最小二乘法是拟合近似曲线，使偏差平方和最小。一般设：
φ(x)=a0φ0(x)+a1φ1(x)+?+a nφn(x)
m
S(a0,a1,……,a n)=∑[a0φ0(x)+a1φ1(x)+?+a nφn(x)?y i]2
i=1
求S的最小值，即
?S
=0?a0(φk,φ0)+a1(φk,φ1)+?+a n(φk,φn)=(φk,f)
k
其中k=0,1…n
对于代数多项式拟合，可以得到：
源程序zxecf.m：
function[a]=zxecf(x,y,n)
% input: x is independent variable matrix
% y is dependent variable matrix
% n is the order of the fitting
% output:a is the polynomial coefficient matrix of fitting
m=length(x);
C=zeros(m,n+1);
for i=1:n+1
C(:,i)=x'.^(i-1);
end
A=C'*C;
b=C'*y';
a=gauss(A,b); %Gaussian elimination
for i=1:n+1
t(i)=a(n-i+2);
end
a=t;
上面代码中用到的高斯消去法的代码如下：
源程序gauss.m：
function [x]=gauss(a,b)
% input the A,B of AX=B
% output the X of AX=B
n=length(a);
ab = [a b]; % augmented matrix
for k=1:n-1
[~,index] = max(abs(ab(k:n,k))); %gain the location of max number
index=index+k-1; % The biggest column main entry number of rows
if index ~= k
temp = ab(index,:);
ab(index,:) = ab(k,:);
ab(k,:) = temp;
end%exchange the column main entry to current row
for j=k+1:n
factor=ab(j,k)/ab(k,k);
ab(j,k:n+1) = ab(j,k:n+1) - factor.*ab(k,k:n+1);
end
end%elimination process
x(n)=ab(n,n+1)/ab(n,n);
for i=n-1:-1:1
sum=ab(i,n+1);
for j=i+1:n
sum=sum-ab(i,j)*x(j);
end
x(i)=sum/ab(i,i);
end
end
以2012年数据为例，输入：
先试进行四阶拟合，运行源程序a=zxecf(x,y,4)，可得：
常数项不为0，可见拟合效果不是很好，继续增大拟合阶数至7阶时，可得：
即y=?6.0662x7+20.2785x6?26.4926x5+17.7751x4?6.6249x3+ 1.9313x2+0.1988x
此时拟合效果比较好。拟合效果如下图：
x=[0 0.05 0.1 0.2 0.4 0.6 0.8 0.9 1]
y=[0 0.0139 0.0342 0.0851 0.2214 0.4034 0.6464 0.8077 1]
a =
0.6832 -0.9742 1.0406 0.2511 -0.0004
a =
-6.0662 20.2785 -26.4926 17.7752 -6.6249 1.9313 0.1988 -0.0000
对所有年份进行七阶拟合，结果如下表所示：
表5.最小二乘法拟合结果
2.龙贝格积分
龙贝格积分法是通过用若干个积分近似值来推算更精确的新的近似值的方法。外推公式为：
R(j,k)=R(j,k?1)+R(j,k?1)?R(j?1,k?1)
4k?1
其中第一列：
T0(0)=b?a
2
[f(a)+f(b)]
T0(l)=1
T0(l?1)+
b?a
l
∑f[a+(2i?1)
b?a
l
]
2l?1
i=1
,l=1,2,3…
外推到k=4即可，因为进一步外推，外推值趋向于0，此时意义不大，所以可以终止进一步的外推。
代码romberg.m：
function [i,r,gi]=romberg(a,n)
% input:a is the poly coefficient matrix
% n is the order of integral
% output:i--integral value
% r is romberg matrix
% gi is the Gini coefficient
f=@(x1)polyval(a,x1);
x=0;
y=1;
h=y-x; % h is the intal steplength
r=zeros(n+1,n+1); % initialize the romberg matrix
r(1,1)=h*(f(x)+f(y))/2;
for j=2:n+1; % calculate the first line
h=h/2;
r(j,1)=r(j-1,1)/2+h*sum(f(x+h*(1:2:2^(j-1))));
end
for k=2:n+1 % calculate the romberg matrix
r(k:n+1,k)=r(k:n+1,k-1)+(r(k:n+1,k-1)-r(k-1:n,k-1))/(4^(k-1)-1);
end
i=vpa(r(n+1,n+1),5);
gi=(0.5-i)/0.5;
将上面得到的2012年的系数矩阵输入，可以得到结果：
i =
0.34238
r =
0.5001 0 0 0 0
0.3909 0.3545 0 0 0
0.3568 0.3455 0.3449 0 0
0.3462 0.3426 0.3424 0.3424 0
0.3433 0.3424 0.3424 0.3424 0.3424 gi =
0.31524714285717436723643913865089
可以得到2012年我国城镇居民的基尼系数为0.3152。
同样经计算可得：
表6.基尼系数计算结果
3.simpson 积分法
Simpson 公式法是采用曲线连接各采样点，而不是像梯形公式一样使用直线连接，然后再将面积求和。由于曲线更逼近于原函数曲线，所以通常情况下，Simpson 公式法的误差要小于梯形公式法。计算公式为：
S =?
3
(f (a )+4∑f (x 2i+1)m?1i=0
+2∑f (x 2i )m?1
i=1
+f(b))
本例不使用最小二乘法所拟合出的公式，直接采用插值积分的方法进行计算。
源代码simpson.m ：
以2012年数据为例，输入：
function [i,gi]=simpson(x,y)
% composite simpon
% input:x,y--the independent and dependent varible matrix % output:i--integral value % gi--the Gini coefficient
n=length(x);
h=(x(n)-x(1))/(n-1);
i=h/3*(y(1)+y(n)+4*sum(y(2:2:n-1))+2*sum(y(3:2:n-2)));
i=vpa(i,5);% 5 effective figures gi=(0.5-i)/0.5;
x=[0 0.05 0.1 0.2 0.4 0.6 0.8 0.9 1]
y=[0 0.0139 0.0342 0.0851 0.2214 0.4034 0.6464 0.8077 1]
可得：
i =
0.33518
gi =
0.32963333333327682339586317539215
即2012年城镇基尼系数为0.3296。
同样经计算可得：
表7.基尼系数计算结果2
四、分析及总结：
前文提到，大陆统计局局长马建堂称，按照国际新的统计口径，大陆居民收入的基尼系数，2003年是0.479，2004年是0.473，2005年为0.485，2006年为0.487，2007年为0.484，2008年为0.491，2009年为0.490，2010年为0.481，2011年为0.477，到2012年的数据是0.474，为2005年以来最低水平。而本文中所计算出的基尼系数却基本都位于0.28-0.40之间，这两组数据差距比较大。分析可得其中有以下原因造成：
1.本文中所采用的数据是中华人民共和国国家统计局的国家数据网中的中国城镇居民的收入指标，并不能够代表全国的人民收入。经查找，数据网上面并没有将农村和城市一并计算的统计表格，所以无法得知国民总平均收入和差距分层等等。农村相对比城镇居民收入应该会更低一些，如果将其计入可能会造成与本计算结果有很大偏差。
2.对居民分层比较粗略，每10%甚至20%分为一等，造成数据点量比较少，所以也可能造成一定的误差。
3.由于我们不知道数据网站中的抽样是如何进行的，所以也并不知道其数据的真实性有几成。比如看2012年城镇居民困难户(5%)的收入，平均为
7,520.90元/年。也就是说，他们平均每个人每个月的收入都有接近650元。这很不符合常理，困难户基本只靠最低生活保障金来生活，即使是在2013年上调过后，一般城市的最低生活保障金也才过400元，更不要说家里面有无收入的儿童。而对于某些富人，其资产收入可能不计其数，而这也是统计局所不能统计到的，所以数据的真实性也有待考据。
4.马建堂本人也说过：“关于大家关心的基尼系数的计算和发布，2012年我曾经说过，中国居民的基尼系数的计算和发布需要城乡住户调查从城乡分开的、城乡收入概念不一致的调查制度，走向中国统一的城乡可比的住户调查制度。也就是说，基尼系数是反映中国居民的收入差异情况，要计算它，就需要中国居民的收入是多少，分等份的收入是多少。过去城乡分开的住户调查，大家也注意到了，只有分城乡的农村居民人均纯收入和城镇居民人均可支配收入，没有中国居民的可支配收入，没有可比的同样指标的城乡居民的收入。”此次国家统计局公布的数据引来舆论较多质疑。各界普遍认为，其数据存在水
分，不足以说明当前收入分配差距过大的社会现实。甚至有极端观点批评称：“童话都不敢这么写”。这可能也造成了本文中的数据差距。
虽然计算出的结果有很大的差距，但是这也给了我很多启发。我国现阶段收入差距还是不小，对此，国家应该改变现行税制在调节收入分配方面的制度缺陷，完善税收调节体系，使税收调节分配的功能在居民收入、存量财产、投资收益等各个环节得到有效发挥、运用综合调控手段，加强对高收入阶层的税收调控、把“富民优先”作为经济发展新阶段以及解决基尼系数拉大问题的重大经济政策，对低收入者实施积极的税收扶持政策以及完善配套措施，加大对非常态高收入阶层收入的监管。
当然，在解决贫富悬殊、化解基尼系数“越警”方面，税收的作用毕竟是有限的，必须和政府其他宏观经济政策一起共同发挥作用，才能更好地解决中国收入分配差距扩大的问题，从而促进中国经济社会健康和谐发展。
参考文献：
【1】中华人民共和国国际统计局，国家数据网https://www.docsj.com/doc/431687675.html,/index
【2】百度百科，基尼系数词条https://www.docsj.com/doc/431687675.html,/view/186.htm
【3】计算方法[M]. 浙江大学出版社, 1989.
【4】程永宏. 改革以来全国总体基尼系数的演变及其城乡分解Ξ[J]. 中国社会科学, 2007 (4).
【5】王金南, 逯元堂, 周劲松. 基于GDP 的中国资源环境基尼系数分析[J][J]. 中国环境科学, 2006, 26(1): 111-115.
数值计算方法大作业
目录 第一章非线性方程求根 (3) 1.1迭代法 (3) 1.2牛顿法 (4) 1.3弦截法 (5) 1.4二分法 (6) 第二章插值 (7) 2.1线性插值 (7) 2.2二次插值 (8) 2.3拉格朗日插值 (9) 2.4分段线性插值 (10) 2.5分段二次插值 (11) 第三章数值积分 (13) 3.1复化矩形积分法 (13) 3.2复化梯形积分法 (14) 3.3辛普森积分法 (15) 3.4变步长梯形积分法 (16) 第四章线性方程组数值法 (17) 4.1约当消去法 (17) 4.2高斯消去法 (18) 4.3三角分解法 (20)
4.4雅可比迭代法 (21) 4.5高斯—赛德尔迭代法 (23) 第五章常积分方程数值法 (25) 5.1显示欧拉公式法 (25) 5.2欧拉公式预测校正法 (26) 5.3改进欧拉公式法 (27) 5.4四阶龙格—库塔法 (28)
数值计算方法 第一章非线性方程求根 1.1迭代法 程序代码： Private Sub Command1_Click() x0 = Val(InputBox("请输入初始值x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = (Exp(2 * x0) - x0) / 5 If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求f(x)=e2x-6x=0在x=0.5附近的根（ep=10-10）
1.2牛顿法 程序代码： Private Sub Command1_Click() b = Val(InputBox("请输入被开方数x0")) ep = Val(InputBox(请输入误差限ep)) f = 0 While f = 0 X1 = x0 - (x0 ^ 2 - b) / (2 * b) If Abs(X1 - x0) < ep Then Print X1 f = 1 Else x0 = X1 End If Wend End Sub 例：求56的值。（ep=10-10）
数值分析上机作业
数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：
课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量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 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。
数值分析上机作业
昆明理工大学工科研究生《数值分析》上机实验 学院：材料科学与工程学院 专业：材料物理与化学 学号：2011230024 姓名: 郑录 任课教师：胡杰
P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误！未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? （1）用MA TLAB函数“eig”求矩阵全部特征值。 （2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。解：MA TLAB程序如下： 求矩阵A的特征值： clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果： 求矩阵B的特征值： clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果：
求矩阵错误！未找到引用源。的特征值： clear; 错误！未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。) 输出结果： （2）A= 10 7877565861097 5 9 10 第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到： 第二部：[Q1,R1]=qr(A1);A2=R1*Q1
matlab与数值分析作业
数值分析作业(1) 1：思考题（判断是否正确并阐述理由） （a）一个问题的病态性如何，与求解它的算法有关系。 (b)无论问题是否病态，好的算法都会得到它好的近似解。 (c)计算中使用更高的精度，可以改善问题的病态性。 (d)用一个稳定的算法计算一个良态问题，一定会得到他好的近似解。 (e)浮点数在整个数轴上是均匀分布。 (f)浮点数的加法满足结合律。 (g)浮点数的加法满足交换律。 (h)浮点数构成有效集合。 (i)用一个收敛的算法计算一个良态问题，一定得到它好的近似解。√2: 解释下面Matlab程序的输出结果 t=0.1; n=1:10; e=n/10-n*t 3：对二次代数方程的求解问题 20 ++= ax bx c 有两种等价的一元二次方程求解公式
2224b x a c x b ac -±==- 对 a=1，b=-100000000，c=1，应采用哪种算法？ 4：函数sin x 的幂级数展开为： 357 sin 3!5!7! x x x x x =-+-+ 利用该公式的Matlab 程序为 function y=powersin(x) % powersin. Power series for sin(x) % powersin(x) tries to compute sin(x)from a power series s=0; t=x; n=1; while s+t~=s; s=s+t; t=-x^2/((n+1)*(n+2))*t n=n+2; end
（a ） 解释上述程序的终止准则； （b ） 对于x=/2π、x=11/2π、x =21/2π，计算的精度是多少？分别需 要计算多少项？ 5:指数函数的幂级数展开 2312!3!x x x e x =+++ + 根据该展开式，编写Matlab 程序计算指数函数的值，并分析计算结果（重点分析0x <的计算结果）。
数值计算方法作业
数值计算方法作业 姓名：李琦 学号：062410124 求 013=--x x 在x=1.5附近的一个根。 一．牛顿下山法： #include
二．加权法 #include
三．单点弦法： #include
(完整版)数值计算方法上机实习题答案
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;
东南大学数值分析上机作业汇总
东南大学数值分析上机作业 汇总 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII
数值分析上机报告 院系： 学号： 姓名：
目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) （1）最大δ值文件 (4) （2）验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序： (10) 2、数据输入及计算结果 (12)
作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序： 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 0.737577 3、两种方法有效位数对比
数值分析作业
第二章 1. 题目：运用MATLAB编程实现牛顿迭代 2. 实验操作 1、打开MATLAB程序软件。 2、在MATLAB中编辑如下的M程序。 function [p1,err,k,y]=newton(f,df,p0,delta,max) %f 是要求根的方程(f(x)=0); %df 是f(x)的导数； %p0是所给初值,位于x*附近； %delta是给定允许误差； %max是迭代的最大次数； %p1是newton法求得的方程的近似解； %err是p0的误差估计； %k是迭代次数； p0 for k=1:max p1=p0-feval('f',p0)/feval('df',p0); err=abs(p1-p0); p0=p1; k p1 err y=feval('f',p1) if (err
p0 = 1.2000 k =1 p1=1.1030 err=0.0970 y=0.0329 k= 2 p1=1.0524 err=0.0507 y=0.0084 k =3 p1=1.0264 err=0.0260 y=0.0021 k =4 p1=1.0133 err=0.0131 y=5.2963e-004 k =5 p1=1.0066 err=0.0066 y=1.3270e-004 k =6 p1=1.0033 err=0.0033 y=3.3211e-005 k =7 p1=1.0017 err=0.0017 y=8.3074e-006 k =8 p1=1.0008 err=8.3157e-004 y = 2.0774e-006 k =9 p1=1.0004 err=4.1596e-004 y =5.1943e-007 k=10 p1=1.0002 err=2.0802e-004 y= 1.2987e-007 k=11 p1=1.0001 err=1.0402e-004 y =3.2468e-008 k=12 p1=1.0001 err=5.2014e-005 y=8.1170e-009 k=13 p1=1.0000 err=2.6008e-005 y= 2.0293e-009 k=14 p1=1.0000 err=1.3004e-005 y=5.0732e-010 k=15 p1 =1.0000 err=6.5020e-006 y=1.2683e-010 k=16 p1 =1.0000 err=3.2510e-006 y=3.1708e-011 k=17 p1 =1.0000 err=1.6255e-006 y =7.9272e-012 k=18 p1 =1.0000 err =8.1279e-007 y= 1.9820e-012 ans = 1.0000 结果说明：经过18次迭代得到精确解为1，误差为8.1279e-007。
数值计算方法第4次作业
第四章 问题一 一、问题综述 在离地球表面高度为y处的重力加速度如下： 计算高度y=55000m处的重力加速度值。 二、问题分析 以高度y作为自变量，重力加速度的值为因变量。得到以下信息： f(0)=9.8100; f(30000)=9.7487; f(60000)=9.6879; f(90000)=9.6278; f(120000)=9.5682; 本题要求的就是f(55000)的值。 以下将采用课堂中学到的Lagrange插值多项式法、Newton插值多项式法、分段低次插值法和样条插值法求解该问题。 三、问题解决 1. lagrange插值多项式法 对某个多项式函数，已知有给定的k+ 1个取值点： 其中对应着自变量的位置，而对应着函数在这个位置的取值。 假设任意两个不同的x j都互不相同，那么应用拉格朗日插值公式所得到的拉格朗日插值多项式为：
其中每个为拉格朗日基本多项式（或称插值基函数），其表达式为： 拉格朗日基本多项式的特点是在上取值为1，在其它的点上取值为0。 源程序lagrange.m function [c,f]=lagrange(x,y,a) % 输入：x是自变量的矩阵；y是因变量的矩阵；a是要计算的值的自变量； % 输出：c是插值多项式系数矩阵；f是所求自变量对应的因变量； m=length(x); l=zeros(m,m); % l是权矩阵 f=0; for i=1:m v=1; for j=1:m if i~=j v=conv(v,poly(x(j)))/(x(i)-x(j)); % v是l_i(x)的系数矩阵 end end l(i,:)=v; % l矩阵的每一行都是x从高次到低次的系数矩阵 end c=vpa(y*l,10); % 对应阶次的系数相加，乘以y，显示10位有效数字 for k=1:m f=f+c(k)*a^(m-k); end 输入矩阵 x=[0 30000 60000 90000 120000] y=[9.81 9.7487 9.6879 9.6278 9.5682] a=55000 再运行源函数，可得： c = [ -2.057613169e-23, 4.938271605e-18, -3.703703702e-14, -0.000002046111111, 9.81] f = 9.6979851723251649906109417384537
东南大学-数值分析上机题作业-MATLAB版
2015.1.9 上机作业题报告 JONMMX 2000
1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果
1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序
数值分析Matlab作业
数值分析编程作业
2012年12月 第二章 14.考虑梯形电阻电路的设计，电路如下： 电路中的各个电流{i1,i2,…,i8}须满足下列线性方程组： 12 123 234 345 456 567 678 78 22/ 2520 2520 2520 2520 2520 2520 250 i i V R i i i i i i i i i i i i i i i i i i i i -= -+-= -+-= -+-= -+-= -+-= -+-= -+= 这是一个三对角方程组。设V=220V，R=27Ω，运用追赶法，求各段电路的电流量。Matlab程序如下： function chase () %追赶法求梯形电路中各段的电流量 a=input('请输入下主对角线向量a='); b=input('请输入主对角线向量b='); c=input('请输入上主对角线向量c='); d=input('请输入右端向量d='); n=input('请输入系数矩阵维数n='); u(1)=b(1); for i=2:n l(i)=a(i)/u(i-1); u(i)=b(i)-c(i-1)*l(i); end y(1)=d(1); for i=2:n y(i)=d(i)-l(i)*y(i-1); end x(n)=y(n)/u(n); i=n-1; while i>0 x(i)=(y(i)-c(i)*x(i+1))/u(i); i=i-1; end x 输入如下:
请输入下主对角线向量a=[0,-2,-2,-2,-2,-2,-2,-2]; 请输入主对角线向量b=[2,5,5,5,5,5,5,5]; 请输入上主对角线向量c=[-2,-2,-2,-2,-2,-2,-2,0]; 请输入方程组右端向量d=[220/27,0,0,0,0,0,0,0]; 请输入系数矩阵阶数n=8 运行结果如下： x = 8.1478 4.0737 2.0365 1.0175 0.5073 0.2506 0.1194 0.0477 第三章 14.试分别用(1)Jacobi 迭代法；(2)Gauss-Seidel 迭代法解线性方程组 1234510123412191232721735143231211743511512x x x x x ?????? ??????---????????????=--?????? --?????? ??????---?????? 迭代初始向量 (0)(0,0,0,0,0)T x =。 （1）雅可比迭代法程序如下： function jacobi() %Jacobi 迭代法 a=input('请输入系数矩阵a='); b=input('请输入右端向量b='); x0=input('请输入初始向量x0='); n=input('请输入系数矩阵阶数n='); er=input('请输入允许误差er='); N=input('请输入最大迭代次数N='); for i=1:n for j=1:n if i==j d(i,j)=a(i,j); else d(i,j)=0; end end end m=eye(5)-d\a; %迭代矩阵 g=d\b; x=m*x0+g; k=1; while k<=N %进行迭代 for i=1:5 if max(abs(x(i)-x0(i))) >er x=m*x+g; k=k+1;
数值计算大作业
数值计算大作业 题目一、非线性方程求根 1.题目 假设人口随时间和当时人口数目成比例连续增长，在此假设下人口在短期内的增长建立数学模型。 （1）如果令()N t 表示在t 时刻的人口数目，β 表示固定的人口出生率，则人口数目满足微分方程() ()dN t N t dt β=，此方程的解为0()=t N t N e β； （2）如果允许移民移入且速率为恒定的v ，则微分方程变成() ()dN t N t v dt β=+， 此方程的解为 0()=+ (1) t t v N t N e e βββ -； 假设某地区初始有1000000人，在第一年有435000人移入，又假设在第一年年底该地区人口数量1564000人，试通过下面的方程确定人口出生率β，精确到 410-；且通过这个数值来预测第二年年末的人口数，假设移民速度v 保持不变。 435000 1564000=1000000(1) e e βββ + - 2.数学原理 采用牛顿迭代法，牛顿迭代法的数学原理是，对于方程0)(=x f ，如果) (x f 是线性函数，则它的求根是很容易的，牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程0)(=x f 逐步归结为某种线性方程来求解。 设已知方程0)(=x f 有近似根k x (假定0)(≠'x f )，将函数)(x f 在点k x 进行泰勒展开，有 . ))(()()(???+-'+≈k k k x x x f x f x f 于是方程0)(=x f 可近似地表示为 ))(()(=-'+k k x x x f x f 这是个线性方程，记其根为1k x +，则1k x +的计算公式为
数值分析上机作业1-1
数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p （E1-1） 显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简 单的）。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19 x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容： 为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数 u ＝roots （a ） 其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根，而函数 b=poly(v) 的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（E1-2）中的扰动项改成18 x ε或其他形式，实验中又有怎样的现象出现？
数值分析作业答案
数值分析作业答案 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange插值基底。 （3）用Newton基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为: ， 所以： 所以f(x)的二次插值多项式为： （2）用Lagrange插值基底 Lagrange插值多项式为： 所以f(x)的二次插值多项式为： (3) 用Newton基底: 均差表如下： xk f(xk) 一阶均差二阶均差 1 0 -1 -3 3/2 2 4 7/ 3 5/6 Newton插值多项式为： 所以f(x)的二次插值多项式为： 由以上计算可知，三种方法得到的多项式是相同的。 6、在上给出的等距节点函数表，若用二次插值求ex的近似值，要使截断误差不超过10-6，问使用函数表的步长h应取多少？ 解：以xi-1,xi,xi+1为插值节点多项式的截断误差，则有 式中 令得 插值点个数
是奇数，故实际可采用的函数值表步长 8、，求及。 解：由均差的性质可知，均差与导数有如下关系： 所以有： 15、证明两点三次Hermite插值余项是 并由此求出分段三次Hermite插值的误差限。 证明：利用[xk，xk+1]上两点三次Hermite插值条件 知有二重零点xk和k+1。设 确定函数k(x): 当或xk+1时k(x)取任何有限值均可； 当时，，构造关于变量t的函数 显然有 在[xk,x][x,xk+1]上对g(x)使用Rolle定理，存在及使得 在，，上对使用Rolle定理，存在，和使得 再依次对和使用Rolle定理，知至少存在使得 而，将代入，得到 推导过程表明依赖于及x 综合以上过程有： 确定误差限： 记为f(x)在[a,b]上基于等距节点的分段三次Hermite插值函数。在区间[xk，xk+1]上有 而最值 进而得误差估计： 16、求一个次数不高于4次的多项式，使它满足，，。
数值计算方法大作业
题目利用数值计算方法求取基尼系数 姓名与学号 指导教师 年级与专业 所在学院
一、问题综述： 基尼系数（Gini coefficient），是20世纪初意大利学者科拉多·吉尼根据劳伦茨曲线所定义的判断收入分配公平程度的指标。是比例数值，在0和1之间。基尼指数（Gini index）是指基尼系数乘100倍作百分比表示。在民众收入中，如基尼系数最大为“1”，最小等于“0”。前者表示居民之间的收入分配绝对不平均（即所有收入都集中在一个人手里，其余的国民没有收入），而后者则表示居民之间的收入分配绝对平均，即人与人之间收入绝对平等，但这两种情况只出现在理论上；因此，基尼系数的实际数值只能介于0～1之间，基尼系数越小收入分配越平均，基尼系数越大收入分配越不平均。 设右图中的 实际收入分配曲线 （红线）和收入分 配绝对平等线（绿 线）之间的面积为 A，和收入分配绝 对不平等线（蓝 线）之间的面积为 B，则表示收入与 人口之间的比例的基尼系数为 A A+B 。 如果A为零，即基尼系数为0，表示收入分配完全平等（红线和绿线重叠）；如果B为零，则系数为1，收入分配绝对不平等（红线和蓝线重叠）。该系数可在0和1之间取任何值。实际上，一般国家的收入分配，既不是完全平等，也不是完全不平等，而是在两者之间，劳伦茨曲线为一条凸向横轴的曲线。收入分配越趋向平等，劳伦茨曲线的弧度越小（斜度越倾向45度），基尼系数也越小；反之，收入分配越趋向不平等，劳伦茨曲线的弧度越大，那么基尼系数也越大。
基尼系数的调节需要国家通过财政政策进行国民收入的二次分配，例如对民众的财政公共服务支出和税收等，从而让收入均等化，令基尼系数缩小。 基尼系数由于给出了反映居民之间贫富差异程度的数量界线，可以较客观、直观地反映和监测居民之间的贫富差距，预报、预警和防止居民之间出现贫富两极分化。因此得到世界各国的广泛认同和普遍采用。 联合国有关组织规定： ●若低于0.2表示收入平均； ●0.2-0.3表示相对平均； ●0.3-0.4表示相对合理； ●0.4-0.5表示收入差距大； ●0.6以上表示收入差距悬殊。 2013年1月18日，中国国家统计局一次性公布了自2003年以来十年的全国基尼系数。大陆统计局局长马建堂称，按照国际新的统计口径，大陆居民收入的基尼系数，2003年是0.479，2004年是0.473，2005年为0.485，2006年为0.487，2007年为0.484，2008年为0.491，2009年为0.490，2010年为 0.481，2011年为0.477，到2012年的数据是0.474，为2005年以来最低水平，而自2008年起，基尼系数也在逐年下降。而此前西南财大调查数据显示，中国的2012年的基尼系数为0.61，但无论是民间统计的数据还是官方统计的数据，结果都遭到学术界质疑，仍具有争议性。 本文将根据网络上国家统计局的数据，利用上面给出的公式来计算我国从2002年以来的城镇居民基尼系数，并将计算出的数据与现有数据进行比较。 全球基尼系数
数值计算方法计算习题
1.已知ln( 2.0)=0.6931;ln(2.2)=0.7885,ln(2.3)=0.8329, 试用线性插值和抛物插值计算.ln2.1的值并估计误差（牛顿插值和拉格朗日插值） 2.已知函数y=sinx 的数表如下,分别用前插和后插公式计算sin0.57891的值，并估算误差。 i x 0.4 0.5 0.6 0.7 )(i x f 0.38942 0.47943 0.56464 0.64422 3. 已知 i x -2 -1 0 1 2 )(i x f 4 2 1 3 5 求)(x f 的二次拟合曲线)(2x p ，并求)0(f '的近似值。 4. 数值积分公式形如 ?'+'++=≈1 )1()0()1()0()()(f D f C Bf Af x S dx x xf 试确定参数D C B A ,,,使公式代数 精度尽量高；（2）设]1,0[)(4 C x f ∈，推导余项公式?-=1 ) ()()(x S dx x xf x R ，并估计 误差。 5. 已知数值积分公式为： )] ()0([)]()0([2)(''20 h f f h h f f h dx x f h -++≈? λ，试确定积分公式中的参数 λ，使其代数精确度尽量高，并指出其代数精确度的次数。
6. 用复化Simpson 公式计算积分 ()? =1 0sin dx x x I 的近似值，要求误差限为 5105.0-?。 7. 已知012113,,4 2 4 x x x ===，给出以这3个点为求积节点在[]0.1上的插值型求积公 式。 8. 给出 900,cos ≤≤x x 的函数表，步长 )60/1(1='=h ，若函数具有5位有 效数字，研究用线性插值求x cos 近似值时的总误差界。 9. 求一个次数不高于4次的多项式)(x P ，使它满足0)0()0(='=P P ， 1)1()1(='=P P ，1)2(=P 。 10. 单原子波函数的形式为bx ae y -=，试按照最小二乘法决定参数a 和b ，已 知数据如下： X 0 1 2 4 y 2.010 1.210 0.740 0.450 11. 分别用梯形公式和辛普森公式计算下列积分：? +1 02 4dx x x 。并估算误差。 12. 用矩阵的克劳特和克利特尔三角分解法求解方程组：??????? ??=??????? ????????? ??7173530103421101002014321x x x x
数值计算方法实习作业模板小
2.1函数图形与极限 2.1.1 实验目的 1.熟悉Mathematica 基本绘图语句。 2.掌握函数极限的有关操作命令。 3.学会利用Mathematica 软件对函数进行分析研究。 4.熟悉Mathematica 二元函数绘图语句。 2.1.2 实验内容 【基本语句】 1.Plot[f[x],{x,xmin,xmax},选项]; 功能: 画出函数f[x] 从min 到max 间的图形; 2.Plot[{f1[x],f2[x],...},{x,xmin,xmax},选项]; 功能: 在同一坐标系下画出函数f1,f2,...的图形。 3. ParametricPlot[{fx,fy},{t,tmin,tmax}]; 功能: 画出参数方程fx=x(t),fy=y(t)的图形; ParametricPlot[{{f1x,f1y},{f2x,f2y}},{t,tmin,tmax}]; 功能：在同一坐标系下画出用参数方程表示的两幅函数图形。 【备注】fx,fy 的给出方式: ⑴fx=x(t) , fy=y(t) ⑵fx=x ,fy=f(x)与fx=f(x) ,fy=x 构成反函数的图形关系 ⑶r=r(t) , fx=r(t)Cos(t) , fy=r(t)Sin(t) 4. Show[tu1,tu2]功能：将tu1及tu2两幅函数图形重叠在一起，将两个函数图形一起显示。 5. Plot3D[f[x,y],{x,x0,x1},{y,y0,y1}] 功能：作出函数f[x,y]在区域[x0,x1]×[y0,y1]上的图形； ParametricPlot3D[{x[u,v],y[u,v],z[u,v]},{u,u0,u1},{v,v0,v1}] 功能：作出参数方程表示的曲面。 6. Limit[f[x],x->x0] 功能：求函数f[x]在x0处的极限。 7. Limit[f[x],x->x0,Direction->+1] 功能：求函数f[x]在x0处的左极限。 8. Limit[f[x],x->x0,Direction->-1] 功能：求函数f[x]在x0处的右极限。 9. Limit[f[x],x->Infinity] 功能：求函数f[x]在 x->无穷时的极限。 10. Limit[f[x],x->-Infinity] 功能：求函数f[x]在 x->负无穷时的极限。 【实验2.1】画出以下函数的图形。 （1）x y ln = 其中]10,1.0[∈x 。 （2）)6 cos(,sin 21π+==x y x y ，其中]6,4[-∈x 。 （3）14233221,,,--====x y x y x y x y ，其中]4,4[-∈x 。 【实验2.2】画出以下函数的图形。 （1）? ??==t y t x sin 其中],0[π∈t 。 （2）???==?? ???==???==t y t x t y t x t y t x 和2 2,其中]2,2[-∈t 。 （3）???===t r y t r x t r sin cos 2cos 9且其中]4,4[ππ-∈t 。 （3）Mathematica 语句： 【实验2.4】利用图形显示命令作出下列函数的图形：
数值分析上机第四次作业
数值分析上机第四次作业 实验项目：共轭梯度法求解对称正定的线性方程组 实验内容：用共轭梯度法求解下面方程组 (1) 123421003131020141100155x x x x -?????? ? ? ?--- ? ? ?= ? ? ?-- ? ? ?-???? ?? 迭代20次或满足()(1) 1110k k x x --∞-<时停止计算。 (2) Ax b =,A 是1000阶的Hilbert 矩阵或如下的三对角矩阵， A[i,i]=4，A[i,i-1]=A[i-1,i]=-1，i=2,3,..,n b[1]=3, b[n]=3, b[i]=2，i=2,3,…,n-1 迭代10000次或满足()()710k k r b Ax -=-≤时停止计算。 （3）*考虑模型问题，方程为 222222(),(,)(0,1)(0,1)(0,)1,(1,),01(,0)1,(,1),01 xy y x u u x y e x y D x y u y u y e y u x u x e x ??+=+∈=???==≤≤==≤≤ 用正方形网格离散化，若取1/,10h N N ==，得到100n =的线性方程组，并用共轭梯度法（CG 法）求解，并对解作图。 实验要求：迭代初值可以取01(,1,...,)ij u i j N ==，计算到32||||10k r -≤停止．本 题有精确解(,)xy u x y e =，这里k u 表示以k ij u 为分量的向量， u 表示在相应点(,)i j 上取值作为分量的向量． 实验一： （1） 编制函数子程序CGmethod 。 function [x,k]=CGmethod(A,b) n=length(A);x=zeros(n,1);r=b-A*x;rho=r'*r; k=0; while rho>10^(-12) & k<20 k=k+1; if k==1 p=r; else beta=rho/rho1; p=r+beta*p; end