数值代数上机作业第七次
二维Ising模型的相变现象
谭昌汇00401146
北京大学数学科学学院科学与工程计算系，100871
完成于2007年6月19日
摘要(Abstract)
本报告主要讨论了二维Ising模型中的内能和比热容随温度变化的相变现象。通过对内能和比热容图像的突变情况，可以确定相变的临界温度βc。
为解决二维Ising模型，本文使用了Monte-Carlo方法，利用Metropolis算法生成一组随机向量，从而对二维Ising模型进行求解。
关键词：二维Ising模型，相变现象，临界温度，Monte-Carlo方法
1介绍(Introduction)
在当代的物理学中，相变问题是一个很前沿的课题，它在很多实际的问题中都起到了重要的作用。而Ising模型是描述铁磁相变的最具代表性的模型。
对于二维的Ising模型，在理论分析上有重要的成果，即可以求得严格解。但尽管如此，并不能给出物理图像。因此，使用计算机处理Ising模型成为了一种十分重要的手段。Monte-Carlo方法就是其中的一种。它可以求出模型的解，并且还可以给出具体的图像与涨落。下面的所有讨论都是围绕使用Monte-Carlo方法求解二维Ising模型的相变问题而展开的。
2相变问题与Ising模型
(The problem of phase change and Ising model)
2.1相变问题(The problem of phase change)
这里所谓的相变，就是指内能、比热容等参数随温度变化产生的突变。这些参数的突变反映了铁磁系统丧失磁性的相变过程。解决相变问题的关键就是将不同温度下的内能和比热容的图像展示出来，从而找出临界温度。
2.2Ising模型初探(The?rst approach to Ising model)
在统计物理中，人们用晶格模型研究磁性材料的相变，每个格点代表一个电子，它有两种自旋的方式，分别记为1和-1。那么，整个系统的内能和比热容可以有以下公式导出。
U M=1
M
σ
exp{?βH(σ)}
Z M
H(σ)
1
M
H(σ)
1
C M=β2
M
{
σ
H2(σ)
exp{?βH(σ)}
Z M
?[
σ
H(σ)
exp{?βH(σ)}
Z M
]}=
β2
M
[ H2(σ) ? H(σ) 2]
其中，H(σ)=?J
σiσj是能量函数，满足Gibbs分布。具体来说，就是相邻项
的乘积之和。比如对于下图，S1状态只需要和与其相邻的S2至S5状态乘积加和
即可。
2.3遍历性(Ergodicity)
在前面的公式中提到了 H(σ) 是能量状态的空间平均，是比较难于计算的。然而，我们可以使用Monte-Carlo方法的思想，进行下面的估计。
H(σ) ≈1
N
N
i=1
H(σ(i)).
其中，右边的序列{σ(i)}N
i=1是马氏链。左式为空间平均，右式为时间平均，这种
空间平均等于时间平均的结果称之为遍历性。
有遍历性的保证，即可用时间序列估计空间平均，而由Monte-Carlo方法思
想，只要适当选取{σ(i)}N
i=1，就可以对时间平均进行估计，最终得到上面的公式
结论。
2.4Metropolis算法(Metropolis algorithm)
剩下的工作就是去生成一组合适的{σ(i)}N
i=1，我们可以使用Metropolis算法来
达到这个目的。这里仅给出算法，不进行进一步的说明，更多的信息可参见参考文献[1]。
[算法]Metropolis算法
步1：根据相应的预选规则由σ(n)产生预选态σ .
步2：计算A=min{1,exp(β?H)}，其中?H=H(σ )?H(σ).
步3：生成均匀分布的随机数r U[0,1].
步4：如果r≤A，则σn+1=σ ；否则σn+1=σn，转步1.
所谓的预选态，一般有两种，一是等可能预选，即在非σ(n)状态中随机选择一个状态。另一种是单点变化预选，变化σn的一个格点取值，总共M（格点数）种不同情况，随机取一个即可。
使用Metropolis算法可直接调用metropolis.m文件即可。其中Hsigma参数指的是H(σ)，是为了加快运算速度而加入的（通过后面的分析可以发现这样的优化很有必要），读入时直接写作f uncH(sigma)即可。
2
2.5小结(Brief summary)
总结一下，讨论二维Ising 模型的相变现象，首先使用Metropolis 算法生成时齐马氏链，然后利用遍历性和Monte-Carlo 算法估计空间平均，最后根据公式求得不同温度下的内能和比热容，作出图形，找到突跃点，即相变的临界点。
理论已经完备，下面的工作就是进行实际的数值实验。
3
数值实验(Numerical experiment)
3.1
生成能量函数(Creating power function)
前面已经大致提到了能量函数选取Gibbs 分布，即紧邻格点的求和。对于边界条件，需要做一些必要的讨论。
我们选取的边界条件是周期条件，即对M ×M 格点来说，第一列（行）和最后一列（行）的自旋方向是对应相等的。要生成这样的分布，只需要先生成一个任意的(M ?1)×(M ?1)的分布，再扩张得到一个满足条件的Gibbs 分布。具体可见下表。
σ11···σ1(M ?1).........σ(M ?1)1···σ(M ?1)(M ?1)
=?
σ11···σ1(M ?1)σ11............σ(M ?1)1···σ(M ?1)(M ?1)σ(M ?1)1σ11···σ1(M ?1)σ11 对(M ?1)×(M ?1)的分布使用Metropolis 算法生成一个状态序列，然后
对序列中的每一个状态进行上面的扩充，最终可以生成符合周期条件要求的{σ(i )}N i =1，可以按规则求解能量函数。
关于能量函数的生成，可以使用程序包中的funcH.m 函数进行求解。函数输入格式σ（(M ?1)×(M ?1)），输出对应的能量值H 。具体的使用方法请参见程序声明。
3.2最初的实验结果(Results at the ?rst time)
一切就绪，先来看一组实验结果。运行Ising.m 函数（具体参数选取方法见函数声明）。选取M =10,N =100，并采取等可能预选。β取0至1间0.1为间隔的一组数，分别求出U 和C ，对β作图。
0.2
0.4
0.6
0.8
1
?2.5
?2
?1.5
?1
?0.5
0.2
0.4
0.6
0.8
1
?1?0.8
?0.6?0.4
?0.20
0.20.4
0.60.8
1这样的图形完全让人摸不着头脑。再来看看使用单点变化预选得到的图形。
3
00.20.40.60.81
?2.5
?2
?1.5
?1
?0.5
000.20.40.60.81
0.1
0.2
0.3
0.4
0.5
0.6
0.7
显然，使用单点变化预选的结果要好于等可能预选。
为什么等可能预选的结论会如此糟糕呢？这就要注意到?H 的值了。由于等可能预选会有可能改变很多格点的值，导致?H 可能会相当大，这使A =min {1,exp (β?H )}≈0，从而极大可能会保留原来的σn 作为σn +1。这样会生成一串相同的σ列，结论当然糟糕。
为了证明这个论点。取M =3，此时由于格点数不多，因此，A 不会太小，等可能预选可能可以得到不错的结果。下面就是一组得到的图形组。
00.20.40.60.81
?3
?2.5
?2
?1.5
?1
?0.5
00.20.40.60.81
0.1
0.2
0.3
0.4
0.5
0.6
0.7
反观单点变化预选，由于只有一个格点的情况发生了改变，无论M 多大，?H 的值都不会太大。因此，在后面的讨论中都使用单点变化预选法。
草草的看一下图，可以发现C 在0.2到0.4之间有较大的跳跃，对U 来讲，也有部分斜率较大，这有可能成为我们所希望得到的临界点。当然，这个图并不能完全说明问题。为了找出最能说明结论的图形，我们需要对实验结果进行优化。
3.3实验结果的优化(Optimize the Results)
要想让实验结果更加出色，一个最原始但却是十分有效的方法就是多次尝试。因为Monte-Carlo 方法和Metropolis 算法都涉及到了随机数，因此每次运行的结果都不一样，多次尝试可以从中选取比较令人满意的值，达到优化实验的结果。
当然，不停尝试多少显得有些不上档次，下面讨论一些提高总体图形表现力手段。它们确实有用，但在此之前，需要先声明的是，由于随机的存在，不可能出现完美的理想结果，有少数不大支持结论的项是很正常的。
首先，提高N 值可以使得Monte-Carlo 算法的近似更好，从而起到优化的作用。举例来看，如果将前面假设中N 取为1000的话，有以下的图形，可以说是个更加令人满意的结果。
4
00.20.40.60.81?2.5
?2
?1.5
?1
?0.5
000.20.40.60.81
0.1
0.2
0.30.4
0.50.60.7
0.8除了N 之外，还有其它的办法优化实验结果。我们的实验都是从全部格点为一的初始态开始的，在前面的若干步中，由于每次仅能改换一个数字，因此随机性不强。为避免这样的问题，可以将前N 1项去掉，不进入时间序列估计的范畴。另外，由于相邻两项有很高的关联度，因此，我们可以每隔N 2项取一个状态记录，而不是每项都记录。这样也许会有好处。
下面就以N 1=800,N 2=10得到的图像。
0.2
0.4
0.6
0.8
1
?2.5
?2
?1.5
?1
?0.5
0.5
0.2
0.4
0.6
0.8
1
00.1
0.20.30.4
0.50.6
0.70.80.9
1这个图像的效果又有了提升，当然，计算速度更慢，因为原来只需要迭代到1000项的σ就可以了，现在需要800+1000×10=18000项才能得到最终的解。看来，一个更好的结果需要付出更大的代价。
我们无法否认的是，无论N 1和N 2如何取值，各项之间还是有关系的，所以，要想减小随机偏离，还有一个很基本的方法，就是多个结果求均值。当然，这又是更费时间的一种策略，但的确是有效果的。同时兼顾时间的复杂性，我们可以在重复次数和N 1,N 2的取值上进行一些取舍和兼顾，生成一组混合策略。比如，取N =1000,N 1=500,N 2=5,times =10，可以得到以下图像。
00.20.40.60.81?2.5
?2
?1.5
?1
?0.5
0.5
00.20.40.60.81
0.1
0.2
0.30.4
0.5
0.60.7
0.8
0.95
这组图像的中的U 的突变更好的显示了出来。
总结一下，要得到一幅好的图像，我们可以通过牺牲时间来达到。采取的措施包括增大时齐列的长度，除去开头的若干项，隔项选取抽样点，多个结果求均值等。需要注意的是，无论用什么样的方法，由于具有随机性，每次得到的结果都不一样，因此，多次尝试并选取较好的图是必要的。
3.4实验结果的加细(Thinning the Results)
通过上面的实验结果，我们可以大致估计β临界值大概在0.4的周围。但由于我们仅仅计算了十个不同温度下的内能和比热容的值，无法更进一步的估计临界值。要想更好的估计，需要对试验结果进行加细，如β取0至1之间以0.02为间隔的一组数。选取策略N =1000,N 1=800,N 2=10,times =1进行计算，得到下面的图形。
00.20.40.60.81
?2.5
?2
?1.5
?1
?0.5
0.5
00.20.40.60.81
0.2
0.4
0.6
0.8
1
1.2
1.4
加细后的图像变得震荡起来，这是很正常的，因为增加了温度点，这点可能偏离理想值，从而形成一个尖角。由于温度点的增多，要想通过不断尝试使每个点都恰好取到理想的地方是相当困难的，因此，我们只有一方面牺牲更多的时间，一方面不去追求完美，而只是观察大致的走向，得到最终的结果。取N =3000,N 1=2000,N 2=10,times =1得到下面的图形。
00.20.40.60.81
?2.5
?2
?1.5
?1
?0.5
0.5
00.20.40.60.81
0.2
0.4
0.6
0.8
1
1.2
1.4
这组策略运行时间相当长，达到数分钟，但得到的结果已经相当令人满意了，毕竟，一共有51个随机产生的点啊。
3.5格点边长的增加(Increase the length of the points)
前面讨论了M =10的情形，现在看一看扩大M 后的结果。方法完全一样，但是值得注意的是，格点增多后，相应的N ,N 1,N 2都应该增大以达到相同的效果。当然，这会耗费更长的时间。因此，这里只做出一个较为粗糙的结果，仅仅作为一种参考。取β间隔0.1，M =100,N =1000,N 1=800,N 2=5,times =1。
6
0.2
0.4
0.6
0.8
1
?2.2
?2?1.8
?1.6?1.4?1.2?1
?0.8?0.60
0.2
0.4
0.6
0.8
1
01
2
34567
89104临界点温度(The temperature in the critical point)
4.1
求解策略(Strategy of the solvation)
由于每张图都不大一样，因此，通过图像只能确定一个大概的临界点温度，比如，β=0.4就是一个大概的点。但由于加细后的震荡变得严重，我们很难说服自己将一幅图的结论作出最终的结果。因此，我们需要大量的计算结果。
我们通过下面的策略来求临界点的温度。选取β从0到1以0.02为间隔。选取参数策略N =500,N 1=500,N 2=3,times =1。选取得到的51个C 值中的最大点对应的β作为临界点。连续运行50次，取平均值即为临界点温度。
4.2数值结果(Numerical results)
尽管设定的参数值都不太大，但由于要循环50次，因此速度仍然比较慢。特别是当M 较大时，运行速度相当慢（M =200是运行了近3小时）。因此，我们没有足够的时间对所有的M 进行实验，仅选取M =10,20,50,100,200几个比较特别的数进行了实验。
下面的组图是使用上面的策略得到的50个β估计值的分布（蓝线）以及平均的β值（红线）。这可以由solvebeta.m 函数导出。由于速度相当的慢，在实际操作中，根据大致的结果估计，β仅选取了0至0.5中的0.02为间隔的总共26个值。
M =10
M =20
7
M=50M=100
M=200
得到的β的估计值总结为下表。
M102050100200
βM0.39080.33440.20040.15520.1356
从上面的图像和图表结论可以发现，随着M的增大，临界点β的值会发生变化，并且趋势是逐渐下降。如果认定啊β=(k B T)?1的话，临界点温度是逐渐上升的。
另外，我们还可以发现，M越大，利用上述策略求得的值的振荡越小，结果越具有可信度（方差很小）。当M→∞时，会有一个极限的值β∞，这个值是多少尚不知道，需要进行更多步骤的数值实验才能得到。当然，根据上面的数据，可以发现β100和β200的值已经相当接近了，可见下降的趋势在减弱，因此可以大胆估计β∞距离β200不会太远。参考文献[4]中称β∞=1/8，是有能够让人接受的。
β∞反映了全空间范围下的铁磁相变临界点的位置，对于很多物理方面的研究是相当有帮助的。
5总结(Summary)
用二维Ising模型模拟铁磁系统的相变是相当成功的。我们得到了很好的图形结果，并且对于不同边长的系统分别求出了值得信服的临界点值，并进行了
8
对比。如果有更多的时间，还可以进行更多次的数值实验，从而得到关于临界点的更多有用的性质。
9
A参考文献(Bibliography)
[1]张平文、李铁军，《数值分析》，2007年1月，北京大学出版社
[2]T E XGuru，《Manual of L A T E X2 》，网络中文版
[3]程卫国、冯峰、王雪梅、刘艺，《Matlab5.3精要、编程及高级应用》（网络
版），2000年4月，机械工程出版社
[4]孙前芳，”二维依辛模型相变附近临界指数β的计算机模拟计算”，《科学技
术与工程》，2005.17，p1235-1237
[5]林旭升，”二维依辛模型相变临界点温度的模拟计算”，《大学物理》，2000.5，p13-15
B软件支持(Software support)
B.1数值计算工具(Numerical calculation tools)
本次报告主要使用的数值计算工具为Matlab。通过调用m文件可以进行相
关的计算。函数使用方法可见函数的声明。所有函数在电子版中附上。
B.2作图工具(Mapping tools)
本次报告的图像是通过Matlab绘制的。图像的最终格式为*.eps，由于此
格式不支持中文导出，故图像中的文字均用了英文批注。部分图形生成后
用Illustrator进行了少量修改。所有图像文件在电子版中附上。
B.3报告排版工具(Typesetting tools)
本次报告使用了L A T E X进行排版，最终电子版以只读型pdf格式上交，同时附
有源文件*.tex供参考和修改。
C感谢(Thanks)
本次报告能够顺利完成，非常感谢蔡啸天同学提供他的电脑，夜以继日的
运行程序。特别是在求作β100、β200等大型问题时，往往需要超过1小时的连续
运行。感谢他和他的电脑作出的努力，让我得到了一个比较好的结论。
10
Contents
1介绍(Introduction)1 2相变问题与Ising模型
(The problem of phase change and Ising model)1
2.1相变问题(The problem of phase change) (1)
2.2Ising模型初探(The?rst approach to Ising model) (1)
2.3遍历性(Ergodicity) (2)
2.4Metropolis算法(Metropolis algorithm) (2)
2.5小结(Brief summary) (3)
3数值实验(Numerical experiment)3
3.1生成能量函数(Creating power function) (3)
3.2最初的实验结果(Results at the?rst time) (3)
3.3实验结果的优化(Optimize the Results) (4)
3.4实验结果的加细(Thinning the Results) (6)
3.5格点边长的增加(Increase the length of the points) (6)
4临界点温度(The temperature in the critical point)7
4.1求解策略(Strategy of the solvation) (7)
4.2数值结果(Numerical results) (7)
5总结(Summary)8 A参考文献(Bibliography)10 B软件支持(Software support)10
B.1数值计算工具(Numerical calculation tools) (10)
B.2作图工具(Mapping tools) (10)
B.3报告排版工具(Typesetting tools) (10)
C感谢(Thanks)10
11
1.谈谈你对该算法的理解：（简单谈一下你是如何理解该算法的？） 先对84阶矩阵进行LU分解，通过Gauss消元法 对下三角形方程组利用前代法解出y，在对上三角方程组 用回代法解出x…. 2.实验内容 function [ L,U ] = LUfac( A ) for k=1:n-1 A(k+1:n,k)=A(k+1:n,k)/A(k,k); A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-A(k+1:n,k)*A(k,k+1:n); end L=tril(A,0); for i=1:n L(i,i)=1; end U=triu(A,0); End //进行LU分解 function [ b ] = TSL( L,b ) n=size(L,1); for j=1:n-1 b(j)=b(j)/L(j,j); b(j+1:n)=b(j+1:n)-b(j)*L(j+1:n,j); end b(n)=b(n)/L(n,n); end //利用前代法解出y function [ b ] = TSU( U,b ) n=size(U,1); for j=n:-1:2 b(j)=b(j)/U(j,j); b(1:j-1)=b(1:j-1)-b(j)*U(1:j-1,j); end b(1)=b(1)/U(1,1); end //利用回代法解出x
主函数程序 A=eye(84); A=6*A; for i=2:84 A(i,i-1)=8; A(i-1,i)=1; End //生成84阶的矩阵A b=ones(84,1); b=b*15; b(1)=7; b(84)=14; [L,U]=LUfac(A);//调用函数LUfac对矩阵A进行分解 y=TSL(L,b);//调用函数TSL求解y x=TSU(U,y); //调用函数TSU求解X 经过matlab…有 x’ ans = 1.0e+008 * Columns 1 through 7 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 8 through 14 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 15 through 21 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 22 through 28 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 29 through 35 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 36 through 42 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 43 through 49 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 Columns 50 through 56 0.0000 0.0000 0.0000 0.0000 0.0000 -0.0000 0.0000 Columns 57 through 63
学号：姓名：成绩： 程序设计基础（C）第七章数组习题 一、选择题 1.下列叙述中错误的是（）。 A）对于double 类型数组，不可以直接用数组名对数组进行整体输入或输出 B）数组名代表的是数组所占存储区的首地址，其值不可改变 C）在程序执行中，数组元素的下标超出所定义的下标范围时，系统将给出“下标越界”的出错信息D）可以通过赋初值的方式确定数组元素的个数 2.下列关于字符串的叙述中正确的是（）。 A）C 语言中有字符串类型的常量和变量 B）两个字符串中的字符个数相同时才能进行字符串大小的比较 C）可以用关系运算符对字符串的大小进行比较 D）空串一定比空格打头的字符串小 3.当用户要求输入的字符串中含有空格时，应使用的输入函数是（）。 A）scanf( ) B）getchar( ) C）gets( ) D）getc( ) 4.若有定义语句：int a[3][6];，按在内存中的存放顺序，a 数组的第10 个元素是（）。 A）a[0][4] B）a[1][3] C）a[0][3] D）a[1][4] 5.已有定义：char a[ ]="xyz",b[ ]={'x', 'y', 'z'};，下列叙述中正确的是（）。 A）数组a 和b 的长度相同 B）a 数组长度小于b 数组长度 C）a 数组长度大于b 数组长度 D）上述说法都不对 6.下列程序的输出结果是（）。 main( ) { char a[7]="a0\0a0\ 0"; int i,j; i=sizeof(a); j=strlen(a); printf("%d %d\n",i,j); } A）2 2 B）7 6 C）7 2 D）6 2 7.下列能正确定义一维数组的选项是（）。 A）int a[5]={0,1,2,3,4,5}; B）char a[ ]={0,1,2,3,4,5}; C）char a={'A', 'B', 'C'}; D）int a[5]="0123"; 8.有以下程序 #include main() { int a[]={2,3,5,4},i; for(i=0;i<4;i++) switch(i%2) { case 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 四、运行结果： 实验二：
数值代数课程设计实验报告 姓名： 班级： 学号： 实验日期： 一、实验名称 代数的数值解法 二、实验环境 MATLAB7.0 实验一、平方根法与改进平方根法 一、实验要求： 用熟悉的计算机语言将不选主元和列主元Gasuss 消元法编写成通用的子程序，然后用编写的程序求解下列方程组 ?????????? ????????????=????????????????????????? ? ? ? ? ? ? ?? ?????? ???? ?--?1415151515768 168 168 168 1681612321 n n n n n x x x x x x 用所编的程序分别求解40、84、120阶方程组的解。 二、算法描述及实验步骤 GAuss 程序如下： （1）求A 的三角分解：LU A =; （2）求解b y =L 得y ； （3）求解y x =U 得x ； 列主元Gasuss 消元法程序如下： 1求A 的列主元分解：LU PA =; 2求解b y P L =得y ； 3求解y x =U 得x ；
三、调试过程及实验结果： %----------------方程系数---------------- >> A1=Sanduijiaozhen(8,6,1,40); >> A2=Sanduijiaozhen(8,6,1,84); >> A3=Sanduijiaozhen(8,6,1,120); >> b1(1)=7;b2(1)=7;b3(1)=7; >> for i=2:39 b1(i)=15; end >> b1(40)=14; >> for i=2:83 b2(i)=15; end >> b2(40)=14; >> for i=2:119 b1(i)=15; end >> b3(120)=14; %----------------方程解---------------- >> x11=GAuss(A1,b1') >> x12=GAuss Zhu(A1,b1') >> x21=GAuss(A2,b2') >> x22=GAuss Zhu(A3,b3') >> x31=GAuss(A3,b3') >> x32=GAuss Zhu(A3,b3') 运行结果：（n=40） GAuss消元法的解即为 x11 = 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 1.0000 列主元GAuss消元法的解即为x12 =
第二次作业及参考答案 1如何设计实验条件使欲了解的电极基本过程（如电化学反应过程）成为电极系统的受控过程？ 答：设计实验条件使欲了解的电极基本过程成为电极系统的受控过程，需要了解该电极过程在电极总过程中的具体位置。例如对于简单电极过程，电极过程等效电路为： 要使电化学反应过程（等效电路元件为Rct）成为电化学测量过程中电极过程的受控步骤，即要使得电极过程的等效电路简化为 则应该设计如下实验条件： （1）采用鲁金毛细管、高导电率的支持电解质或断电流法、恒电位仪正反馈法等措施，以减小或补偿掉电解液欧姆电阻，电极的电子导体选用高导电率材料作电子导电物质，以减小或忽略掉电子导体的欧姆电阻； （2）电化学测量时采用小幅值外加激励信号，当外加激励作用于电极一段时间、双电层充（放）电结束但浓差极化还没出现时进行测量，以忽略双电层充放电过程和浓差极化的影响。 （3）当电化学反应物可溶时，可采用旋转圆盘电极、在适宜的高转速下对电极进行测量，以忽略浓差极化对电极过程研究的干扰。 2什么是支持电解质？作用是什么？ 答：支持电解质：电导率强、浓度大、在电化学测量过程中承担溶液电迁移任务而不参与电化学反应的物质。可以使溶液的离子强度和电导率在测量过程中保持稳定，获得重现性良好的可靠数据。作用：（1）增强溶液导电性，减小溶液欧姆压降；（2）承担溶液电迁移任务，消除反应物或产物的电迁移传质；（3）支持电解质浓度大，离子迁移数大，溶液离子强度主要取决于支持电解质，可以忽略测量过程中因反应物或产物离子消耗引起的离子强度的变化，电极反应各物种扩散系数可近似视为常数；（4）有利于紧密双电层电容的构建，减小分散层电势（1电势）的影响；（5）加入支持电解质后溶液密度大，可以忽略因电活性物质浓度梯度引起的溶液密度差，从而减小或忽略界面附近的对流传质。 3 工作电极分类？ 答：按电极是否作为反应物参与电极反应，工作电极分为两类：第一类工作电极和第二类工作电极。第一类工作电极可称为活性工作电极，电极既承担电子导电的任务，又作为反应物参与电极反应。第二类工作电极又称为惰性工作电极，
中国海洋大学本科生课程大纲 课程属性：公共基础/通识教育/学科基础/专业知识/工作技能，课程性质：必修、选修 一、课程介绍 1.课程描述： 数学实验是由于计算机技术和科学计算软件的迅猛发展应运而生的一门较新的数学课程，它改变了数学只靠纸和笔的传统形象，将实验的手段引入到数学的学习和研究中。 本课程为大学二年级数学院的学生开设。它不是讲授新的数学知识，而是让学生利用已有的数学知识去解决一些经简化的实际问题。大多数实验的一般过程是：对于给出的实际问题，建立数学模型、选择适当的数学方法、用科学计算软件MATLAB编程计算、对运算结果进行分析、给出结论。 本课程以MATLAB软件为主要的实验工具，采用以学生动手动脑为主，教师讲授和点评、小组讨论、报告为辅的教学方式。 通过本课程的学习，学生用数学解决实际问题的意识和能力可以得到强化和提高，更切实地体会到数学的用处，增加学习兴趣，提高创造力。 2.设计思路： 本课程旨在训练用数学解决实际问题的能力。实验内容的选取是基于学生具备MATLAB语言的初步编程能力、并学习了数学分析、高等代数、解析几何、运筹学基础（初步）、数学实验基础、常微分方程、数值分析或计算方法、概率论等数学课程的基础之上。课程共分七个基础实验和一个综合实验依次进行。七个基础实验是：MATLAB 基础知识复习、常微分方程（组）、数据建模——插值与拟合、古典密码学、图与网络 - 6 -
优化、动态规划、遗传算法。 基础实验涉及的数学内容较为单一、数学模型和求解方法较简单，是对“用数学”能力的基本训练。 综合实验以三人为一组进行，所涉及到的数学知识范围更广，建模和求解的难度更大。综合实验的题目可以小组自拟或在任课教师拟定的题目中选择。任课教师拟定的题目将于综合实验开始前一周给出。各小组在实验前要上交一份“开题报告”：写出问题的重述、模型建立和求解的思路、可能遇到的主要困难及解决方案。通过认真完成综合实验，“用数学”的能力可以有一个较大的提升。 3.课程与其他课程的关系： 先修课程：高等代数I、高等代数II、空间解析几何、数学分析I、数学分析II、数学实验基础；常微分方程；计算方法（或数值分析、数值代数）； 并行课程：概率论等； 后置课程：数学模型；数学建模实践 二、课程目标 本课程的目标是为大二数学类专业学生提供用数学知识解决实际问题的系统训练。 到课程结束时，学生应能： （1）对简单的实际问题建立数学模型； （2）采用适当的数学方法，用MA TLAB软件求解模型，并根据计算结果对模型进行评价和改进； （3）具备初步的科研写作能力：学会如何将问题、模型、解决思路、求解方法、计算结果和结论简洁、清晰、严谨地呈现； （4）针对难度较高的实际问题通过小组成员的独立思考、相互合作与激励，共同解决。提高沟通交流能力，促进相互学习，加深对有关数学知识的理解，进一步提升用数学知识和MATLAB软件解决实际问题的能力。 三、学习要求 要完成所有的课程任务，学生必须： （1）按时上课,认真听讲，积极参与课堂讨论、随堂练习和测试； - 6 -
1、依次输入五句话，然后将它逆序输出。运行结果如图所示。 提示： 创建一个字符串数组，每句话作为字符串数組的一个元素，然后从该数组的末尾开 始循环输出。 2、某百货商场当日消费积分最高的八名顾客.他们的积分是：18、25、7、36、1 3、2、89、63。编写程序找出最低的积分及它在数组中的原始位置（下标）。 提示： > 创建数组points[]，存储八名顾客的积分。 > 定义变量min存储最低枳分，定义变量index存储最低积分下标。 > 假设第一个元素为最低积分，下标为0。 > 遍历数组，将数组元素和min的值进行比较。 3、从键盘输入10个整数，合法值为1、2或3，不是这三个数则为非法数字。试编程统计每个整数和非法数字的个数。运行结果如图所示。
提示： ?创建数组nums[]，长度为10，存储用户输入的数字。 ?创建数组count[]，长度为4、存储三个合法数字和非法数字的个数。 ?循环输入数字，利用switch判断数字的值、根据不同的值对数组count[]中的不同元素进行累加。 4、假设有一个长度为5的数组，如下所示： int [] array = new int []{1,3,-1,5,2}; 现要创建一个新数组newArray[],要求新数组中的元素与原数组中的元素逆序，并且如果 原数组中的元素值小于0，在新数组中按0存储。试编程输出新数组中的元素，程序运行结果如图： 提示： ?利用循环从原数组最后一个元素（下标为array.length-1）开始处理，如杲该元素的值小于0,利用continue退出本次循环（整型数组中元素默认值为0)。 ?如果该元素值大于0,则将该元素复制到新数组 ?合适的位置。对于原数组下标为i的元素.在新数组中的下标为arrayJength-i-1。 ?处理完成，利用循环输出新数组中的元素 5、有一数列：8, 4, 2, 1, 23, 344, 12。编写程序，让用户输入一个整数，判断该整数在不在该数列中。运行效果如下图所示：
实验报告名称 班级：学号：姓名：成绩： 1实验目的 1）通过对二分法与牛顿迭代法作编程练习与上级运算，进一步体会二分法与牛顿迭代法的不同特点。 2）编写割线迭代法的程序，求非线性迭代法的解，并与牛顿迭代法。 2 实验内容 用牛顿法和割线法求下列方程的根 x^2-e^x=0; x*e^x-1=0; lgx+x-2=0; 3实验步骤 1)根据二分法和牛顿迭代法，割线法的算法编写相应的求根函数； 2)将题中所给参数带入二分法函数，确定大致区间； 3)用牛顿迭代法和割线法分别对方程进行求解； 3 程序设计 牛顿迭代法x0=1.0; N=100; k=0; eps=5e-6; delta=1e-6; while(1) x1=x0-fc1(x0)/fc2(x0); k=k+1; if k>N disp('Newmethod failed')
break end if(abs(x1-x0) 华东理工大学 复变函数与积分变换作业（第5册） 班级____________学号_____________姓名_____________任课教师_____________ 第九次作业 教学内容：5.1孤立奇点 5.2.1 留数的定义 5.2.2极点处留数的计算 1．填空题： （1）函数)1(1)(i e z f z +-= 的全部孤立奇点是 ,.......1,0),24 ( 2ln 2 1±=++k k i ππ （2）0=z 是 z z -sin 1的____三____级极点. （3）2-=z 是 3 2 3) 4(8--z z 的____三_级极点. （4）若()f z 在0z 点解析，0z 是()g z 的本性极点，0z 是()()f z g z ?的_本性_奇点, 是 ()() f z g z 的___本性__奇点. （5）=]0,1cos [Res z z 2 1 2．指出下列函数的奇点及其类型（不考虑∞点），若是极点，指出它的级. (1) 21n n z z +; 解：由,1,01-==+n n z z 得) 1,,1,0()12(-==+n k e z i n k k π 为原式一级极点。 (2) z z ) 1(ln + 解1：10,1 ) 1()1ln(0 1 <<+-= +∑∞ =+z n z z n n n ， ∑∞ =+-= +0 1 ) 1() 1ln(n n n n z z z 无负幂项，故0 =z 为其可去奇点。 解2：1) 1(1lim ) 1ln(lim =+=+→→z z z z z ，故0=z 为可去奇点。 (3)1z z e - 解：由于1z z e -∑∞ =---+---= ==0 11 1 111) 1() 1(1 n n n z z z z e e e e ，所以1=z 为本性奇点。 (4) 3 sin z z ; 解：0=z 为z sin 的一级零点，而0=z 为3z 的三级零点，故0=z 为 3 sin z z 的二级极点。 01sin lim sin lim 0 3 2 ≠==→→z z z z z z z ，故0=z 为 3 sin z z 的二级极点。 (5) 2 1(1) z z e -; 解：因),! 32 1()! 1(12 ++ + =+=-∑ ∞ =z z z n z z e n n z 故0=z 为 2 1(1) z z e -的三级极点，而 ,2,1,2±±==k i k z π均为一级极点。 (6) 2 sin z e z z 解：由于 2 sin z e z z z z e z z z e z z ...) ! 31(....) ! 3(2 2 3 +- = +-= 所以1sin lim 2 =? →z z e z z z ，因此 0z =是一级极点 3 证明：如果0z 是()f z 的(1)m m >级零点，那么0z 是()f z '的1m -级零点. 证明：0z 是()f z 的()1m m >级零点，可设()()()0m f z z z z ?=-， 其中()z ?在0z 点解析，且()00z ?'≠， 科技信息 SCIENCE &TECHNOLOGY INFORMATION 2012年第9期偏微分方程（PDE ）是众多描述物理，化学和生物现象的数学模型的基础，其最新应用已经扩展到经济，金融预测，图像处理等很多领域。要通过PDE 模型研究这些问题，就需要求解PDE 方程，但是绝大多数微分方程特别是偏微分方程，很难得到其解析形式的解。我们希望能够借助于计算机采用数值方法求得偏微分方程的近似解，这就是《偏微分方程数值解》课程的主要内容。 《偏微分方程数值解》是信息与计算科学专业的一门专业课，它与《数值代数》，《数值逼近》一起构成信息与计算科学专业信息与计算方向的核心课程，在专业培养中占有非常重要的地位。随着计算机技术的飞速发展，偏微分方程数值解得到了前所未有的发展和应用，与此同时也暴露了《偏微分方程数值解》课程传统教学中的很多不足之处，这使得该门课程在教学上有很多地方需要调整。 笔者长年教授《偏微分方程数值解》课程，在该门课程的教学改革方面做了一些思索和尝试，主要包括改革教学方法，更新教学模式，加强介绍背景知识，融入数学建模思想，教学与科研相结合，教学与计算软件相结合，增设实验课，改革考核方式等。 1改革教学方法，更新教学模式 由于数学课程大多理论性较强，趣味性较弱，为了激发学生学习兴趣，在教学过程中，我们采用启发式、讨论式等多种教学方法，营造良好的课堂气氛，加强师生之间的交流，引导学生独立思考，强化科学思维的训练。在教学内容方面，不光教授公式推导，定理证明，同时注重算法思想的讲解和程序设计的讲解，同时安排一定课时的习题课，讲解典型习题和对每章进行总结。 由于《偏微分方程数值解》涉及较多的概念、公式和定理，大多数老师仍以传统的课堂教学为主，而少数年轻教师则喜欢用多媒体课件教学。传统的教学方法，虽然受到的批评最多，但也是用得最多，最能让大家普遍接受的一种方法，在算法推导、理论分析等方面，采用传统的板书讲解能更好地引导学生去感受和思考数学逻辑的过程以及创造性的思维过程，加深对数学理论的理解和认识，培养学生的逻辑和思维能力。而在讲述背景知识，算法的应用，算法的程序实现时候最好用多媒体课件进行演示。多媒体课件可以让学生更直观，更全面的理解算法的应用，另外使用多媒体课件还可以节省大段公式的板书时间，图示清楚、准确。但是如果全部使用多媒体课件上课，容易加快教学速度，淡化数学公式的推导以及定理的证明过程，不利于培养学生的数学思维能力。所以，我们认为需要将传统的教学方法和现代的教学手段结合起来，充分发挥各自的优势，在传统教学中穿插多媒体课件，根据教学内容选择合适的教学手段。 2加强知识背景的介绍，融入数学建模思想 《偏微分方程数值解》是理论知识与实际应用之间的桥梁，为学生使用计算机解决科学与工程中的实际问题打下良好的理论基础和应用基础。传统教学以分析，证明，推导为主，重理论，轻应用，缺少偏微分方程产生的实际背景的介绍和应用数值解的方法解决实际问题的实例。因此，我们在教授该课程的时候，注重与数学建模思想相结合，从实际问题出发，建立相应的偏微分方程模型，这样，学生就知道为什么要研究偏微分方程，偏微分方程能解决什么样的实际问题。 例如，我们考虑有衰减的扩散问题：有一个扩散源，某物质从此扩散源向四周扩散，沿x,y,z 三个方向的扩散系数分别为常数，衰减使质量的减少与浓度成正比，扩散前周围空间此物质的浓度为0，估计物质的分布。我们引导学生运用所学过的微积分的思想以及相应的物理知识，对这一问题进行建模，可以得到如下的模型： 鄣u =a 2鄣2 u 鄣x +b 2鄣2 u 鄣y +c 2鄣2 u 鄣z -k 2u 上述方程是常系数线性抛物型方程，它就是有衰减的扩散过程的数学模型。有了这样的铺垫，学生知道了扩散问题的数学模型就是抛物型方程，当然类似的环境污染，疾病流行等与扩散有关的实际问题可以用抛物型方程来描述，很自然的，接下来的问题就是如何求解上面的抛物型方程，学生的学习热情自然就提高了。 3教学与科研相结合 随着计算技术和计算机科学的发展，偏微分方程数值解法的内涵也在不断扩大，我们在讲授《偏微分方程数值解》课程中引进近年来最新的理论和最新的方法，这样可以开阔学生的视野，激发学生的学习情趣，锻炼学生的自学能力。例如我们除了介绍有限差分法，有限元法，有限体积法等经典的具有一般性的方法，还介绍了多重网格法。由于近些年来，人们将辛方法应用于哈密顿常微分方程系统以及推广应用于微分方程的兴趣日益增长，我们也简单介绍了这一主题，并且用这个思想去分析逼近波动方程的交错蛙跳格式。在讲授方法的同时，还注意介绍这些方法的发展历史，设计思想和理论依据，并给出了相当丰富的参考文献，让基础好的同学自己去挖掘感兴趣的问题。承担课题的老师，可以把自己课题中与此课程相关的小问题拿出来供有兴趣的同学琢磨，有助于锻炼学生的科研能力。 4教学与计算软件相结合 由Mathworks 公司推出的MATLAB 软件，现在已经发展成功能强大，适合科学和工程计算的软件，使用MATLAB 编程，语言简洁，数据处理方便，具有强大的数值计算功能和图形展示功能，因此，将MATLAB 融入偏微分方程数值解的教学，更能与时俱进，更有效地提高教学质量。 MATLAB 采用有限元的方法求解各种PDE ，它提供了两种方法解决PDE 问题，一是pdepe 函数，它可以求解一般的PDEs ，具有较大的通用性，但只支持命令行形式的调用。二是PDE 工具箱，可以求解特殊PDE 问题，但有较大的局限性。只能求解二阶PDE 问题，不能求解偏微分方程组。PDE 工具箱支持命令行形式求解，但需要记住大量命令及其调用格式。不过好在它提供了GUI 界面，可以把我们从复杂的编程中解脱出来，还有很好的动画演示功能，尤其适合刚入门的学生。 我们在授课过程中精选与生活，生产密切相关的应用实例，鼓励学生自己动手建立模型，应用数学软件和所学的知识求解模型。例如考虑一个带有矩形孔的金属板上的热传导问题。板的左边保持在100℃，板的右边热量从板向环境空气定常流动，其他边及内孔边界保持绝缘。初始t=t 0时板的温度为0。对于这样的一个实际问题，我们先应用所学的数学分析和数学建模知识，对原问题建立如下偏微分方程模型： 鄣u 鄣t -△u =0，u =100， 鄣u =-1，鄣u =0，u|t=t 0 =0△△△△△△△△△△△△△△ △. 不妨设界顶点坐标为（-0.5，-0.8），（0.5，-0.8），（0.5，0.8），（-0.5，0.8）。内边界顶点坐标为（-0.005，-0.4），（0.05，-0.4），（0.05，0.4），（-0.05，0.4）。对于这样的一个抛物型方程，我们设计其数值计算方法，然后分别用 偏微分方程数值解课程的思索 邹永魁 （吉林大学数学与科学学院吉林 长春 130012） 【摘要】探讨《偏微分方程数值解》课程教学改革的思考与体会，主要包括教学方法和教学模式的改革，加强背景知识的介绍，将科研前沿带入课堂，将MATLAB 融入教学以及考核方式的改革等。 【关键词】偏微分方程数值解；教学改革；MATLAB ；综合评价体系○高校讲坛○200 数值分析上机实验报告 《数值分析》上机实验报告 1.用Newton 法求方程 X 7-X 4+14=0 在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。 1.1 理论依据： 设函数在有限区间[a ，b]上二阶导数存在，且满足条件 {}α?上的惟一解在区间平方收敛于方程所生的迭代序列 迭代过程由则对任意初始近似值达到的一个中使是其中上不变号 在区间],[0)(3,2,1,0,) (') ()(],,[x |))(),((|,|,)(||)(|.4;0)(.3],[)(.20 )()(.110......b a x f x k x f x f x x x Newton b a b f a f mir b a c x f a b c f x f b a x f b f x f k k k k k k ==- ==∈≤-≠>+ 令 )9.1()9.1(0)8(4233642)(0)16(71127)(0)9.1(,0)1.0(,1428)(3 2 2 5 333647>?''<-=-=''<-=-='<>+-=f f x x x x x f x x x x x f f f x x x f 故以1.9为起点 ?? ?? ? ='- =+9.1)()(01x x f x f x x k k k k 如此一次一次的迭代，逼近x 的真实根。当前后两个的差<=ε时，就认为求出了近似的根。本程序用Newton 法求代数方程(最高次数不大于10)在（a,b ）区间的根。 1.2 C语言程序原代码： #include 本科第三次作业 （第十六周课程讲完之后交，用学院统一的作业纸书写，学院文印室有卖，要求抄题。）1.一个系有多个学生，每个学生只能在一个系注册；一个学生可以选修多门课程，每门课 程可以有许多个学生选修。用一个ER图表示“系”、“学生”、“课程”的数据联系。 2.设对乐曲的采样频率是每秒44 100次，采样值用32位表示。计算录制1小时的音乐需 要的多少存储容量？ 存储容量=44100×32×3600/(8×1024×1024)=605.6M 3.显示器的解像度为1024×768位，每个像素的颜色要用16个位来表示，计算一幅画面 需要多少个字节来存储。 显示器的解像度为1024×768位，每个像素的颜色要用16个位来表示，计算一幅画面需要多少个字节来存储。需要字节数=1024×768×16/8=1572864 4.数据管理技术的发展经历了那几个阶段？ 数据管理技术的经历了3个阶段，即：人工管理阶段、文件系统阶段和数据库系统阶段。 5.数据库技术的主要特点是什么？ 数据库系统的主要特点是：(1)数据的结构化、(2)数据的共享性、 (3)数据的独立性、(4)数据的完整性、(5)数据的灵活性、(6)数据的安全性。 6.教材202页第4题。描述波形音频和MIDI音乐的区别？ 波形音频和MIDI音乐的区别：与声音波形相比，MIDI数据不是声音而是指令，所以它的数据量要比波形声音少的多。MIDI可在多媒体应用中与其他波形声音配合使用，形成伴奏的效果。MIDI声音尚不能做到在音质上与真正的乐器完全相似，无法模拟出自然界中其他乐曲类声音。 7.教材202页第6题。多媒体数据为什么可以压缩？视频标准有哪些？ 由于多媒体数据为中的相关性很强，并且有大量的冗余信息，当人们采用适当的压缩方法后，可以得到很大的压缩比。常用的视频标准有：AVI、DVAVI、 8.教材202页第7题。简述数据加密和解密的工作原理是什么？ 加密的工作原理是发送对数据进行伪装，即使这些数据被窃取，非法用户得到的也是一对杂乱无章的垃圾数据，不能获得任何信息。解密的工作原理合法用户接收导数据后，通过事先指定的处理方法将这些数据还原为原始数据。 9.教材202页第8题。什么样的计算机程序被称为计算机病毒？计算机病毒具有哪些特 - 1 - 第三章上机习题 用你所熟悉的的计算机语言编制利用QR 分解求解线性方程组和线性最小二乘问题的 通用子程序，并用你编制的子程序完成下面的计算任务： （1）求解第一章上机习题中的三个线性方程组，并将所得的计算结果与前面的结果相比较，说明各方法的优劣； （2）求一个二次多项式+bt+c y=at 2 ，使得在残向量的2范数下最小的意义下拟合表3.2中的数据； （3）在房产估价的线性模型 111122110x a x a x a x y ++++= 中，1121,,,a a a 分别表示税、浴室数目、占地面积、车库数目、房屋数目、居室数目、房龄、建筑类型、户型及壁炉数目，y 代表房屋价格。现根据表3.3和表3.4给出的28组数据，求出模型中参数的最小二乘结果。 （表3.3和表3.4见课本P99-100） 解 分析： （1）计算一个Householder 变换H ： 由于T T vv I ww I H β-=-=2，则计算一个Householder 变换H 等价于计算相应的v 、β。其中)/(2,||||12v v e x x v T =-=β。 在实际计算中， 为避免出现两个相近的数出现的情形，当01>x 时，令2 12221||||) (-x x x x v n +++= ； 为便于储存，将v 规格化为1/v v v =，相应的，β变为)/(221v v v T =β 为防止溢出现象，用∞||||/x x 代替 （2）QR 分解： 利用Householder 变换逐步将n m A n m ≥?,转化为上三角矩阵A H H H n n 11 -=Λ，则有 ?? ? ???=0R Q A ，其中n H H H Q 21=，:),:1(n R Λ=。 在实际计算中，从n j :1=，若m j <，依次计算)),:((j m j A x =对应的)1()1()~ (+-?+-k m k m j H 即对应的j v ，j β，将)1:2(+-j m v j 储存到),:1(j m j A +，j β储存到)(j d ，迭代结束 后再次计算Q ，有??? ? ?? ??=-~001 j j j H I H ，n H H H Q 21=（m n =时1-21n H H H Q =） （3）求解线性方程组b Ax =或最小二乘问题的步骤为 i 计算A 的QR 分解； ii 计算b Q c T 11=，其中):1(:,1n Q Q = iii 利用回代法求解上三角方程组1c Rx = （4）对第一章第一个线性方程组，由于R 的结果最后一行为零，故使用前代法时不计最后一行，而用运行结果计算84x 。 运算matlab 程序为 1 计算Householder 变换 [v,belta]=house(x) function [v,belta]=house(x) n=length(x); x=x/norm(x,inf); sigma=x(2:n)'*x(2:n); v=zeros(n,1); v(2:n,1)=x(2:n); if sigma==0 belta=0; else alpha=sqrt(x(1)^2+sigma); if x(1)<=0 v(1)=x(1)-alpha; else v(1)=-sigma/(x(1)+alpha); end belta=2*v(1)^2/(sigma+v(1)^2); v=v/v(1,1); end end A.int a[]="string"; B.int a[5]={0,1,2,3,4,5}; C.char a="string"; D.char a[]={0,1,2,3,4,5}; A.1D B.3 C.9 D. A. B.\"c:\\abc.dat\" C."c:\abc.dat" D."c:\\abc.dat" A. B. C. D. A. B. C. D. A. a[2][4] B. a[l C. a[l+l][0] D. a(2)(1) 确定 A. a[0][2*1] B. a[l][3] C. a[4-2][0] D. a[0][4] A. 2 B. 3 C. 4 D. A. 3 5 7 B. 3 6 9 C. 1 5 9 D 1 4 7 A. j*m+i B. i*m+j C. i*m+j-1 D. i*m+j+1 确定 A. B. C. D. A.1D B.3 C.9 D. A. B.\"c:\\abc.dat\" C."c:\abc.dat" D."c:\\abc.dat" A. if(s1>s2) B. if(strcmp(s1,s2)) C. if(strcmp(s2,s1)>O) D. if(strcmp(s1 A. B. C. D. A. int a[3][ ]; B. float a(3 C. double a[1][4]; D. float a(3)(4); A. a[2][4] B. a[l C. a[l+l][0] D. a(2)(1) A. 2 B. 3 C. 4 D. A. 3 5 7 B. 3 6 9 C. 1 5 9 D 1 4 7 A B C D 计算方法实验报告 班级： 学号： 姓名： 成绩： 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.顺序计算 #include9-10次作业答案
偏微分方程数值解课程的思索
数值分析上机实验报告
本科第三次作业答案
数值线性代数二版徐树方高立张平文上机习题第三章实验报告
南昌大学第七章数组答案
计算方法上机实习题大作业(实验报告).
相关主题
文本预览