数值分析大作业

超声波对金属凝固力学性能影响的曲线拟合

摘要：本文采用了超声波功率从0w变化到700w时处理低熔点亚共晶Al-7.3%Si合金，对使用超声波的功率-抗拉强度曲线进行分析。采用了多项式拟合的计算方法对给定的功率和测定的抗拉强度进行拟合。首先采用线性方程对数值进行拟合，计算程序为采用C语言编程计算其参数，计算方法采用多项式进行拟合。然后又采用二次多项式方程对数值进行拟合，计算方法与一次多项式拟合时一样，且采用C语言进行编程。

结果表明: 超声波处理提高了该合金的抗拉强度和延伸率，采用二次多项式进行拟合其误差较小，且能较好的符合功率-抗拉强度曲线之间的变化趋势。

关键词：超声波，力学性能，最小二乘法，多项式拟合，C语言

1 前言

现代科学技术的发展不仅对材料性能要求越来越高,而且对环保要求也日趋严格。材料的环境协调性促使绿色材料加工技术的飞速发展，传统的冶金化学手段细化凝固组织工艺受到了环境理念的质疑和挑战，凝固技术正朝着高效，环保的方向发展。如何能在不污染环境及材料的前提下实现对金属凝固过程和凝固组织的控制是冶金及材料工作者长期追求和奋斗的目标。基于上述发展思路,在凝固过程中施加物理场处理技术成为提高材料性能的重要工艺手段之一。外加物理场处理技术是在金属凝固前或凝固过程中对金属熔体施加物理场，利用金属和物理场相互作用,改善其凝固过程和组织的一种技术，该技术具有环境友好和操作简便等优点。尤其是大功率超声波由于其独特的声学效应对金属凝固过程具有十分显著的影响，因此在材料领域的研究和应用再次成为热点。

其中超声波的功率对材料的力学性能具有重大的现实意义，为工艺参数的选择和工业应用提供了理论依据。但是几个特定的功率并不能很好的反应材料在整个处理过程中的功率-抗拉强度变化。因此就不能研究材料在处理过程中的力学性能变化的情况。通过采用数值分析曲线拟合的知识可以较为简便的近似得出材料在处理过程中的功率-抗拉强度变化曲线，对实际工作中课题研究上影响很大。

2 正文部分

(1) 最小二乘法的基本原理

从整体上考虑近似函数 同所给数据点 (i =0,1,…,m)误差 (i =0,1,…,m)的大小，常用的方法有以下三种：一是误差 (i=0,1,…,m)绝对值的最大值，即误差向量的∞-范数；二是误差绝对值的和，即误差向量r 的1-范数；三是误差平方和的算术平方根，即误差向量r 的2-范数；前两种方法简单、自然，但不便于微分运算，后一种方法相当于考虑 2-范数的平方，因此在曲线拟合中常采用误差平方和来度量误差 (i=0，1，…，m)的整体大小。

在函数的最佳平方逼近中()[]b a C x f ,∈，如果()x f 只在一组离散点集

{}m i x i ???=,1,0,上给出，这就是科学实验中经常可以见到的实验数据

(3) 即

（4）即上式子为法方程组，其是关于的线性方程组，用矩阵表示为

（5）可以证明，方程组(5)的系数矩阵是一个对称正定矩阵，故存在唯一解。从

式(5)中解出(k=0,1,…，n)，从而可得多项式

(6) 可以证明，式（6）中的满足式（1），即为所求的拟合多项式。我

们把称为最小二乘拟合多项式的平方误差，记作

多项式拟合的一般方法可归纳为以下几步：

1)由已知数据画出函数粗略的图形——散点图，确定拟合多项式的次数n；

2)列表计算和；

3)写出正规方程组，求出；

4)写出拟合多项式。

在实际应用中，或；当时所得的拟合多项式就是拉格朗日

或牛顿插值多项式。

(2) 问题的描述

本文首先采用此技术处理了低熔点亚共晶Al-7.3%Si合金, 未处理试样的抗

拉强度为155MPa,200W处理后试样的抗拉强度稍有提高为163MPa,随后400W和600W处理,试样抗拉强度分别为182MPa和ZllMPa"700W超声波处理后,合金的抗拉强度达232MPa,比未处理试样的提高了50%"。试通过用多项式拟合的方法近似的计算合金在处理时的超声波功率-抗拉强度曲线。测试的几个点如表1所示。

表1 合金在超声波功率处理下的不同抗拉强度

为方便计算通过把表1中变换一下为表2

表2 合金在超声波功率处理下的不同抗拉强度

（3）问题的求解过程

首先对数据进行一次多项式拟合即令S 1(x)=a+bx ，φ0(x)=1，φ1(x)=x ，

（φ0，φ0）=

∑=7

1i =1x1+1x1+…1x1=8

（φ1，φ0）=（φ0，φ1）=∑=70

i i x =0+1+2+…+7=28 （φ1，φ1）=∑=7

02i i x =0+1+4+…+49=140

（f ，φ0）=170

?∑=i i y =155+160+163…+232=1470

（f ，φ1）=i i i x y ?∑=7

=0x155+1x160+2x163…+7x232=5597

由法方程组得：

???? ??=???? ?????? ??5597

147014028288

b a 解得 a=146.09 b=10.76

因此其线性方程组为 ()x x S 76.1009.1461+=

然后再对数据进行一次多项式拟合即令S 1(x)=a+bx+cx 2，φ0(x)=1，φ1(x)=x ，φ2(x)=x 2

（φ1，φ2）=（φ2，φ1）=∑=7

03i i x =0+1+8+…+343=784

（φ2，φ2）=∑=7

4i i x =0+1+16+…+22401=4676

（f ，φ2）=27

i i i x y ?∑==0x155+1x160+4x163+…+49x232=29127

（φ1，φ1）=（φ2，φ0）=（φ0，φ2）=140

所以由上式且有法方程组得：

???

?

? ??=????? ??????? ??2912755971470467678414078414028140288

c b a

解得方程组得 a=156.03 b=0.82 c=1.42

所以此二次多项式为 ()2142.182.003.156x x x S ++=

以上结果的计算程序也可以通过C 语言程序计算而来的，计算程序见附录。

（4）问题的计算结果及其分析

由编程程序可得一次多项式的结果如下：

由编程程序可得二次多项式的结果如下：

(5) 图形的对比

由以上计算结果可知：

1)由法方程组计算的一次和二次多项式与编程所得的式子高度一样，其误差可能是由计算机的舍入造成的

2)由计算的均方误差和偏差可知，二次多项式的均方误差和偏差比一次多项式都要小，且图形更加吻合。这就说明用二次多项式拟合超声波功率-抗拉强度关系比一次多项式更好，更能吻合曲线的分布和趋势。

3 参考文献

[1] 李清扬,王能超,易大义.数值分析 [M]. 北京:清华大学出版社,2008,51-96

[2] 谭浩强,C语言程序设计[M].北京:清华大学出版社,1991

[3] 李红,徐长发.数值分析学习辅导.习题解析[M].武汉:华中科技大学出版社，

2006,73-109

完成人签名：

完成时间：2013-12-15

【附录】

1)用C语言编程一次多项式的程序如下：

/*用二维数组实现2X3阶矩阵的三角直接分解(由于系数矩阵是对称正定矩阵，不必换元) a[i][j](i=0,1;j=0,1)用于存放系数矩阵A，a[i][2](i=0,1)用于存放向量b

*/

#include

#include

void main()

{

void juzh(float b[2][3],float f1[],float x1[],int n);

void zhy(float b[2][3],int n);

void jie(float b[2][3]);

void pia(float b[],float f1[]);

float ss(float b[]);

float x[8]={0,1,2,3,4,5,6,7};

float f[8]={155,160,163,171,182,196,211,232};

float s[8],ip[2];

float a[2][3];

int i,j;

for(i=0;i<=1;i++)

juzh(a,f,x,i);

for(i=0;i<=1;i++)/*求比例因子*/

{ip[i]=a[i][0];

for(j=0;j<=1;j++)

{if(fabs(ip[i])

ip[i]=a[i][j];

}

}

for(i=0;i<=1;i++)

{for(j=0;j<=2;j++)

a[i][j]=a[i][j]/ip[i];

}

for(i=0;i<=1;i++)

zhy(a,i);

jie(a);

for(i=0;i<=7;i++)

{s[i]=a[0][2]+a[1][2]*x[i];

}

pia(s,f);

printf("换元三角直接分解系数LU矩阵：解向量x=");

for(i=0;i<=1;i++)

{for(j=0;j<=2;j++)

{if(j%3==0)

printf("\n");

printf("%5.10f ",a[i][j]);

}

}

printf("\n");

printf("拟合曲线表达式为：f(x)=%.6f+%.6fx",a[0][2],a[1][2]);

printf("\n");

printf("偏差为：");

for(i=0;i<=7;i++)

{if(i%4==0) printf("\n");

printf("%.10f ",s[i]);

}

printf("\n");

printf("均方差为：");

printf("%.10f",ss(s));

printf("\n");

}

void juzh(float b[2][3],float f1[],float x1[],int n)/*求内积矩阵（法方程）*/ {int i,j;

float sum=0;

for(i=0;i<=1;i++)

{sum=0;

for(j=0;j<=7;j++)

sum+=pow(x1[j],n)*pow(x1[j],i);

b[n][i]=sum; }

sum=0;

for(i=0;i<=7;i++)

sum+=pow(x1[i],n)*f1[i];

b[n][2]=sum;

}

void zhy(float b[2][3],int n)/*换元三角分解*/

{int i,j;

float sum=0;

if(n==0)

{for(i=n;i<=1;i++)

b[i][n]=b[i][n];

}

else

{for(i=n;i<=1;i++)

{sum=0;

for(j=0;j<=n-1;j++)

sum-=b[i][j]*b[j][n];

b[i][n]=b[i][n]+sum;

}

}

for(i=n+1;i<=1;i++)

b[i][n]=b[i][n]/b[n][n]; sum=0;

if(n==0)

{for(i=n+1;i<=1;i++)

b[n][i]=b[n][i];

}

else

{for(i=n+1;i<=1;i++)

{sum=0;

for(j=0;j<=n-1;j++)

sum-=b[n][j]*b[j][i];

b[n][i]=b[n][i]+sum;

}

}

}

void jie(float b[2][3])

{float sum=0;

int i,j;

for(i=1;i<=1;i++)/*实现Ly=Pb*/ {sum=0;

for(j=0;j<=i-1;j++)

sum-=b[i][j]*b[j][2];

b[i][2]=b[i][2]+sum;

}

sum=0;

b[1][2]=(b[1][2])/(b[1][1]);

for(i=0;i>=0;i--)/*实现Ux=y*/ {sum=0;

for(j=1;j>=i+1;j--)

sum-=b[i][j]*b[j][2];

b[i][2]=b[i][2]+sum;

b[i][2]=b[i][2]/b[i][i];

}

}

void pia(float b[],float f1[])

{int i;

for(i=0;i<=7;i++)

{b[i]=b[i]-f1[i];

}

}

float ss(float b[])

{float sum=0;

int i;

for(i=0;i<=7;i++)

sum=sum+b[i]*b[i];

}

sum=sqrt(sum);

return(sum);

}

2)用C语言编程二次多项式的程序如下：

*用二维数组实现3X4阶矩阵的三角直接分解(由于系数矩阵是对称正定矩阵，可不必换元) a[i][j](i=0,1,2;j=0,1,2)用于存放系数矩阵A，a[i][3](i=0,1,2)用于存放向量b

#include

#include

void main()

{

void juzh(float b[3][4],float f1[],float x1[],int n);

void zhy(float b[3][4],int n);

void jie(float b[3][4]);

void pia(float b[],float f1[]);

float ss(float b[]);

float x[8]={0,1,2,3,4,5,6,7};

float f[8]={155,160,163,171,182,196,211,232};

float s[8],ip[3];

float a[3][4];

int i,j;

for(i=0;i<=2;i++)

juzh(a,f,x,i);

for(i=0;i<=2;i++)/*求比例因子*/

{ip[i]=a[i][0];

for(j=0;j<=2;j++)

{if(fabs(ip[i])

ip[i]=a[i][j];

}

}

for(i=0;i<=2;i++)

{for(j=0;j<=3;j++)

a[i][j]=a[i][j]/ip[i];

}

for(i=0;i<=2;i++)

zhy(a,i);

jie(a);

for(i=0;i<=7;i++)

{s[i]=a[0][3]+a[1][3]*x[i]+a[2][3]*pow(x[i],2);

}

pia(s,f);

printf("换元三角直接分解系数LU矩阵：解向量x="); for(i=0;i<=2;i++)

{for(j=0;j<=3;j++)

{if(j%4==0)

printf("\n");

printf("%5.10f ",a[i][j]);

}

}

printf("\n");

printf("拟合曲线表达式为：f(x)=%.6f+%.6fx+%.6fx^2",a[0][3],a[1][3],a[2][3]); printf("\n");

printf("偏差为：");

for(i=0;i<=7;i++)

{if(i%4==0) printf("\n");

printf("%.10f ",s[i]);

}

printf("\n");

printf("均方差为：");

printf("%.10f",ss(s));

printf("\n");

}

void juzh(float b[3][4],float f1[],float x1[],int n)/*求内积矩阵（法方程）*/

{int i,j;

float sum=0;

for(i=0;i<=2;i++)

{sum=0;

for(j=0;j<=7;j++)

sum+=pow(x1[j],n)*pow(x1[j],i);

b[n][i]=sum; }

sum=0;

for(i=0;i<=7;i++)

sum+=pow(x1[i],n)*f1[i];

b[n][3]=sum;

}

void zhy(float b[3][4],int n)/*换元三角分解*/

{int i,j;

float sum=0;

if(n==0)

{for(i=n;i<=2;i++)

b[i][n]=b[i][n];

}

else

{for(i=n;i<=2;i++)

{sum=0;

for(j=0;j<=n-1;j++)

sum-=b[i][j]*b[j][n];

b[i][n]=b[i][n]+sum;

}

}

j=n;

for(i=n+1;i<=2;i++)

b[i][n]=b[i][n]/b[n][n];

sum=0;

if(n==0)

{for(i=n+1;i<=2;i++)

b[n][i]=b[n][i];

}

else

{for(i=n+1;i<=2;i++)

{sum=0;

for(j=0;j<=n-1;j++)

sum-=b[n][j]*b[j][i];

b[n][i]=b[n][i]+sum;

}

}

}

void jie(float b[3][4])

{float sum=0;

int i,j;

for(i=1;i<=2;i++)/*实现Ly=Pb*/ {sum=0;

for(j=0;j<=i-1;j++)

sum-=b[i][j]*b[j][3];

b[i][3]=b[i][3]+sum;

}

sum=0;

b[2][3]=(b[2][3])/(b[2][2]);

for(i=1;i>=0;i--)/*实现Ux=y*/ {sum=0;

for(j=2;j>=i+1;j--)

sum-=b[i][j]*b[j][3];

b[i][3]=b[i][3]+sum;

b[i][3]=b[i][3]/b[i][i];

}

}

void pia(float b[],float f1[])

{int i;

for(i=0;i<=7;i++)

{b[i]=b[i]-f1[i];

}

}

float ss(float b[])

{float sum=0;

int i;

for(i=0;i<=7;i++)

{

sum=sum+b[i]*b[i];

}

sum=sqrt(sum);

return(sum);

}

数值分析课程学习心得体会

在当今信息时代，数学知识在科学研究、工程技术、人文社会科学以及经济生活等领域中的作用越来越显现重要。可以说，通过数值分析课程的学习，不

仅使我们获取一定的数学算法理论和算法应用能力，而更重要的是培养我们今后的科学素质，开发我们的创新意识，培育我们的创新能力。所以数值分析课程在数学专业和工科有关专业的课程体系中占有十分特殊的地位。为此我对于这一学期的数值分析课程具有很深的体会及感悟。

数值分析课程的内容是研究算法，其基础是数学分析与微分方程（或高等数学），高等代数（或线性代数），要学好数值分析课程，必须对以上的有关内容比较熟悉。数值分析课程内容既有纯数学的高度抽象性、严密的理论性和逻辑性，又有应用的广泛性与实际实验的技术性、近似性。即理论与实践的紧密结合是课程的特点之一。这对于层次不齐的我们来说带来了一定的困难，所以对同一个问题，我们给出的算法可能是不同的，从而应用效果也有一定区别。数值分析课程中的内容都是为了解决实际问题而产生的，所以课程内容具有广泛的物理背景或实际应用价值。其各类算法在许多工程设计、科学计算、社会、经济和生态领域的问题中都有广泛的应用。随着计算技术与科学技术的发展，算法在更多的领域内发挥越来越重要的核心作用。研究算法是为了更好地解决实际问题，因此课程内容对培养我们解决实际问题的能力有很大的作用，对培养我们的创新意识和科学素养有重要的帮助。

在一开始，我就问我们的黄老师，问“为什么要学这个课或内容？学了有什么用”？后来慢慢才发现我们是要真正的理解这些思想，因为这些思想才是真正接近生活的。我们不能满足于能够运用定理而已，我们更要结合实际深刻的理解定理、算法，不断的有意识的无意识的发现并接受定理、算法中蕴含的思想。让这些思想内化为自我知能的一部分，去引领我们的生活。我认为这样的数学才是美的。这样不断有意识无意识的将本已工具化的数学转化为内在的知识，进而真正让数学帮助我们全方位的成长。比如“等价转换”的思想，这里的“等价”不是完全意义上的“等价”，而是要求转换前后转换主体的主要特征没变。插值法的基本思想是抓住已知函数或者已知点的几个主要特征，用另外一个具备了这几个主要特征的简单的函数来代替已知函数或拟合已知数据点。埃尔米特插值就是要求插值多项式在节点上与原函数的函数值，节点上对应的导数值相等；数值积分和数值微分中的主要思想是设法构造某个简单函数 p(x)近似 f(x),然后对p(x)求积（求导）得到f(x)的积分（导数）的近似值。

对于这门课程来说，也给我在精神上带来感悟，主要有一下几点：

第一，兴趣，爱因斯坦说过“兴趣是最好的老师。”我觉得非常有道理，我

曾对于数学很感兴趣，因此在学习过程中不会感到枯燥乏味，一个人一旦对某事物有了浓厚的兴趣，就会主动去求知、去探索、去实践，并在求知、探索、实践中产生愉快的情绪和体验。这样才能主动学习，并且学好到精通。

第二，勤奋，一勤天下无难事。从古到今，有多少名人不是有勤奋而得来成功的。现在的年轻人，一代比一代聪明。要不被别人淘汰，要超越别人，只有靠时间堆出来。每天多学一些，多积累一些。再算大作业的时候，数据是非常繁杂的，这就需要你坚持不懈，不能放弃，一步一步的算下去，一点都不能懈怠。。

第三，悟性，通常人认为指顿悟，慧根，我觉得就是对一个问题不断的思索，将自己的体会和感受融合，获得属于自己的知识。有很多事情、问题，都是可以想明白的。只有不停的想，才能想明白，想透彻。在学习过程中，一些定理和算法你要举一反三的来看，比如在Dolittle分解和Crout分解中，就是可以相互转换的。

数值分析作为大学工科专业的公共课程，其强调和注重它算法的应用过程与方法，。因此我希望老师可以介绍有关理论与算法时，根据课程特点，注重借助一些典型的物理背景或实际问题，从问题出发组织教学内容，以突出课程的应用性来激发学生的学习兴趣，降低具体算法的理论与难度，让我们便于理解算法。对重要的算法还可以根据不同层次学生的情况，布置不同难度或形式的上机实验题目，使每个学生都有自己的实践体会和理解，可以做到因势利导、因材施教。对于考核来说，我觉得在分值分配上是否可以进行一些调整，比如考试成绩占80%，而大作业占20%，因为这样可以避免一些同学不会编程的尴尬。

最后衷心感谢老师的辛勤教育，因为有你们才能让我们接受更多的知识，才能让我们懂得如何更好的去学习。虽然这门课程已经结束，但这不是终点，而是一个崭新的起点。

数值分析大作业-三、四、五、六、七

大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0; end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +??=-= ?-???解：Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while nerrorlim n=n+1; else break ; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数 非线性函数f function y=f(x) y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end

数值分析大作业三 四 五 六 七

大作业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:');

flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag==1 sigma=k*eps; x0=sigma; k=k+1; m=0; flag1=1; while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0;

end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =-= ?-???解： Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1;

数值分析作业答案

数值分析作业答案 插值法 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次的多项式，使它满足，，。

数值计算方法大作业

目录 第一章非线性方程求根 (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）

数值分析大作业

数值分析报大作业 班级：铁道2班 专业：道路与铁道工程 姓名：蔡敦锦 学号：13011260

一、序言 该数值分析大作业是通过C语言程序编程在Microsoft Visual C++ 6.0编程软件上运行实现的。本来是打算用Matlab软间来计算非线性方程的根的。学习Matlab也差不多有一个多月了，感觉自己编程做题应该没什么问题了；但是当自己真心的去编程、运行时才发现有很多错误，花了一天时间修改、调试程序都没能得到自己满意的结果。所以，我选择了自己比较熟悉的C程序语言来编程解决非线性的求值问题，由于本作业是为了比较几种方法求值问题的收敛速度和精度的差异，选择了一个相对常见的非线性函数来反映其差异，程序运行所得结果我个人比较满意。编写C语言，感觉比较上手，程序出现问题也能比较熟练的解决。最终就决定上交一份C程序语言编程的求值程序了！

二、选题 本作业的目的是为了加深对非线性方程求根方法的二分法、简单迭代法、、牛顿迭代法弦截法等的构造过程的理解；能将各种方法的算法描述正确并且能够改编为程序并在计算机上实现程序的正确合理的运行，能得到自己满意的结果，并且能调试修改程序中可能出现的问题和程序功能的增减修改。本次程序是为了比较各种方法在求解同一非线性方程根时，在收敛情况上的差异。 为了达到上面的条件我选择自己比较熟悉的语言—C语言来编程，所选题目为计算方程f(x)=x3-2x-5=0在区间[2,3]内其最后两近似值的差的绝对值小于等于5 ?的根的几种方法的比较。 110- 本文将二分法、牛顿法、简单迭代法、弦截法及加速收敛法这五种方法在同一个程序中以函数调用的方式来实现，比较简洁明了，所得结果能很好的比较，便于分析；发现问题和得出结论。

北航数值分析大作业一

《数值分析B》大作业一 SY1103120 朱舜杰 一．算法设计方案： 1.矩阵A的存储与检索 将带状线性矩阵A[501][501]转存为一个矩阵MatrixC[5][501] . 由于C语言中数组角标都是从0开始的，所以在数组MatrixC[5][501]中检索A的带内元素a ij的方法是： A的带内元素a ij=C中的元素c i-j+2,j 2.求解λ1，λ501，λs ①首先分别使用幂法和反幂法迭代求出矩阵按摸最大和最小的特征值λmax和λmin。λmin即为λs； 如果λmax>0,则λ501=λmax；如果λmax<0,则λ1=λmax。 ②使用带原点平移的幂法（mifa（）函数），令平移量p=λmax，求 出对应的按摸最大的特征值λ，max， 如果λmax>0,则λ1=λ，max+p；如果λmax<0,则λ501=λ，max+p。 3.求解A的与数μk=λ1+k（λ501-λ1）/40的最接近的特征值λik （k=1,2，…，39）。 使用带原点平移的反幂法，令平移量p=μk，即可求出与μk最接近的特征值λik。 4.求解A的（谱范数）条件数cond（A）2和行列式d etA。 ①cond（A）2=|λ1/λn|，其中λ1和λn分别是矩阵A的模最大和 最小特征值。

②矩阵A的行列式可先对矩阵A进行LU分解后，detA等于U所有对角线上元素的乘积。 二．源程序 #include #include #include #include #include #include #include #define E 1.0e-12 /*定义全局变量相对误差限*/ int max2(int a,int b) /*求两个整型数最大值的子程序*/ { if(a>b) return a; else return b; } int min2(int a,int b) /*求两个整型数最小值的子程序*/ { if(a>b) return b; else return a; } int max3(int a,int b,int c) /*求三整型数最大值的子程序*/ { int t; if(a>b) t=a; else t=b; if(t

北航数值分析大作业第二题精解

目标：使用带双步位移的QR 分解法求矩阵10*10[]ij A a =的全部特征值，并对其中的每一个实特征值求相应的特征向量。已知：sin(0.50.2)() 1.5cos( 1.2)(){i j i j ij i j i j a +≠+== (i,j=1,2, (10) 算法： 以上是程序运作的逻辑，其中具体的函数的算法，大部分都是数值分析课本上的逻辑，在这里特别写出矩阵A 的实特征值对应的一个特征向量的求法： ()[]()() []()[]()111111I 00000 i n n n B A I gause i n Q A I u Bu u λλ-?-?-=-?-?? ?-=????→=??????→= ?? ? 选主元的消元 检查知无重特征值 由于=0i A I λ- ，因此在经过选主元的高斯消元以后，i A I λ- 即B 的最后一行必然为零，左上方变 为n-1阶单位矩阵[]()()11I n n -?-，右上方变为n-1阶向量[]()11n Q ?-，然后令n u 1=-，则 ()1,2,,1j j u Q j n ==???-。

这样即求出所有A所有实特征值对应的一个特征向量。 #include #include #include #define N 10 #define E 1.0e-12 #define MAX 10000 //以下是符号函数 double sgn(double a) { double z; if(a>E) z=1; else z=-1; return z; } //以下是矩阵的拟三角分解 void nishangsanjiaodiv(double A[N][N]) { int i,j,k; int m=0; double d,c,h,t; double u[N],p[N],q[N],w[N]; for(i=0;i

上海大学_王培康_数值分析大作业

数值分析大作业（2013年5月） 金洋洋（12721512），机自系 1.下列各数都是经过四舍五入得到的近似值，试分别指出它 们的绝对误差限, 相对误差限和有效数字的位数。 X1 =5.420, x 2 =0.5420, x 3=0.00542, x 4 =6000, x 5=50.610? 解：根据定义：如果*x 的绝对误差限 不超过x 的某个数位的半个单位，则从*x 的首位非零数字到该位都是有效数字。 显然根据四舍五入原则得到的近视值，全部都是有效数字。 因而在这里有：n1=4, n2=4, n3=3, n4=4, n5=1 （n 表示x 有效数字的位数） 对x1：有a1=5, m1=1 （其中a1表示x 的首位非零数字，m1表示x1的整数位数） 所以有绝对误差限 143 11 (1)101022 x ε--≤ ?=? 相对误差限 31() 0.510(1)0.00923%5.4201 r x x x εε-?= == 对x2：有a2=5, m2=0 所以有绝对误差限 044 11 (2)101022 x ε--≤ ?=? 相对误差限 42() 0.510(2)0.00923%0.54202 r x x x εε-?= == 对x3：有a3=5, m3=-2 所以有绝对误差限 235 11 (3)101022 x ε---≤ ?=? 相对误差限 53() 0.510(3)0.0923%0.005423 r x x x εε-?= == 对x4：有a4=0, m4=4 所以有绝对误差限 4411(4)1022 x ε-≤?= 相对误差限 4() 0.5 (4)0.0083%6000 4 r x x x εε= = = 对x5：有a5=6, m5=5 所以有绝对误差限 514 11(5)101022 x ε-≤ ?=? 相对误差限 45() 0.510(5)8.3%600005 r x x x εε?= ==

北航数值分析报告第三次大作业

数值分析第三次大作业 一、算法的设计方案： （一）、总体方案设计： x y当作已知量代入题目给定的非线性方程组，求（1）解非线性方程组。将给定的(,) i i

得与(,)i i x y 相对应的数组t[i][j],u[i][j]。 （2）分片二次代数插值。通过分片二次代数插值运算，得到与数组t[11][21],u[11][21]]对应的数组z[11][21]，得到二元函数z=(,)i i f x y 。 （3）曲面拟合。利用x[i],y[j],z[11][21]建立二维函数表，再根据精度的要求选择适当k 值，并得到曲面拟合的系数矩阵C[r][s]。 （4）观察和(,)i i p x y 的逼近效果。观察逼近效果只需要重复上面（1）和（2）的过程，得到与新的插值节点(,)i i x y 对应的(,)i i f x y ，再与对应的(,)i i p x y 比较即可，这里求解 (,)i i p x y 可以直接使用（3）中的C[r][s]和k 。 （二）具体算法设计： （1）解非线性方程组 牛顿法解方程组()0F x =的解* x ，可采用如下算法： 1）在* x 附近选取(0) x D ∈，给定精度水平0ε>和最大迭代次数M 。 2）对于0,1, k M =执行 ① 计算() ()k F x 和()()k F x '。 ② 求解关于() k x ?的线性方程组 () ()()()()k k k F x x F x '?=- ③ 若() () k k x x ε∞∞ ?≤，则取*()k x x ≈，并停止计算；否则转④。 ④ 计算(1) ()()k k k x x x +=+?。 ⑤ 若k M <，则继续，否则，输出M 次迭代不成功的信息，并停止计算。 （2）分片双二次插值 给定已知数表以及需要插值的节点，进行分片二次插值的算法： 设已知数表中的点为： 00(0,1,,) (0,1,,)i j x x ih i n y y j j m τ=+=???=+=?? ，需要插值的节点为(,)x y 。 1) 根据(,)x y 选择插值节点(,)i j x y ： 若12h x x ≤+ 或12 n h x x ->-，插值节点对应取1i =或1i n =-，

数值分析第二次大作业

《数值分析》计算实习报告 第二题 院系：机械工程及自动化学院 学号： 姓名： 2017年11月

一、题目要求 试求矩阵A =[a ij ]10×10的全部特征值，并对其中的每一个实特征值求相应的特征向量，已知 a ij ={ sin (0.5i +0.2j ) i ≠j 1.52cos (i +1.2j ) i =j (i,j =1,2, (10) 说明： 1.用带双步位移的QR 方法求矩阵特征值，要求迭代的精度水平为ε=10?12。 2.打印以下内容： （1）全部源程序； （2）矩阵A 经过拟上三角化后所得的矩阵A (n?1)； （3）对矩阵A (n?1)实行QR 方法迭代结束后所得的矩阵； （4）矩阵A 的全部特征值λi =(R i ,I i ) (i =1,2,?,10)，其中R i =Re(λi )，I i = Im(λi ) 。若λi 是实数，则令I i =0； （5）A 的相应于实特征值的特征向量。 3.采用e 型数输出实型数，并且至少显示12位有效数字。 二、算法设计思路和方案 1. 将矩阵A 拟上三角化得到矩阵A (n?1) 为了减少计算量，一般先利用Householder 矩阵对矩阵A 作相似变换，把A 化为拟上三角矩阵A (n?1)，然后用QR 方法计算A (n?1)的全部特征值，而A (n?1)的特征值就是A 的特征值。具体算法如下： 记(1)A A =，()r A 的第r 列至第n 列的元素为(r)(1,2, ,;,1,,)ij a i n j r r n ==+。 对于1,2,,2r n =-执行 （1）若() (2,3,,)r ir a i r r n =++全为零，则令(1)()r r A A +=，转(5)；否则转（2）。

北航数值分析大作业第一题幂法与反幂法

《数值分析》计算实习题目 第一题： 1. 算法设计方案 （1）1λ，501λ和s λ的值。 1)首先通过幂法求出按模最大的特征值λt1，然后根据λt1进行原点平移求出另一特征值λt2，比较两值大小，数值小的为所求最小特征值λ1，数值大的为是所求最大特征值λ501。 2)使用反幂法求λs ，其中需要解线性方程组。因为A 为带状线性方程组，此处采用LU 分解法解带状方程组。 （2）与140k λλμλ-5011=+k 最接近的特征值λik 。 通过带有原点平移的反幂法求出与数k μ最接近的特征值 λik 。 （3）2cond(A)和det A 。 1）1=n λλ2cond(A)，其中1λ和n λ分别是按模最大和最小特征值。 2）利用步骤（1）中分解矩阵A 得出的LU 矩阵，L 为单位下三角阵,U 为上三角阵，其中U 矩阵的主对角线元素之积即为det A 。 由于A 的元素零元素较多，为节省储存量，将A 的元素存为6×501的数组中，程序中采用get_an_element()函数来从小数组中取出A 中的元素。 2.全部源程序 #include #include void init_a();//初始化A double get_an_element(int,int);//取A 中的元素函数 double powermethod(double);//原点平移的幂法 double inversepowermethod(double);//原点平移的反幂法 int presolve(double);//三角LU 分解 int solve(double [],double []);//解方程组 int max(int,int); int min(int,int); double (*u)[502]=new double[502][502];//上三角U 数组 double (*l)[502]=new double[502][502];//单位下三角L 数组 double a[6][502];//矩阵A int main() { int i,k; double lambdat1,lambdat2,lambda1,lambda501,lambdas,mu[40],det;

北航数值分析大作业第二题

数值分析第二次大作业 史立峰 SY1505327

一、 方案 （1）利用循环结构将sin(0.50.2)() 1.5cos( 1.2)() {i j i j ij i j i j a +≠+==(i,j=1,2,……,10)进行赋值，得到需要变换的 矩阵A ； （2）然后，对矩阵A 利用Householder 矩阵进行相似变换，把A 化为上三角矩阵A (n-1)。 对A 拟上三角化，得到拟上三角矩阵A (n-1)，具体算法如下： 记A(1)=A ，并记A(r)的第r 列至第n 列的元素为()n r r j n i a r ij ,,1,;,,2,1) ( +==。 对于2,,2,1-=n r 执行 1. 若 ()n r r i a r ir ,,3,2) ( ++=全为零，则令A(r+1) =A(r),转5；否则转2。 2. 计算 () ∑+== n r i r ir r a d 1 2 )( ()( )r r r r r r r r r r d c a d a c ==-=++则取,0sgn ) (,1)(,1若 )(,12r r r r r r a c c h +-= 3. 令 () n T r nr r r r r r r r r R a a c a u ∈-=++) ()(,2)(,1,,,,0,,0 。 4. 计算 r r T r r h u A p /)(= r r r r h u A q /)(= r r T r r h u p t /= r r r r u t q -=ω T r r T r r r r p u u A A --=+ω)()1( 5. 继续。 （3）使用带双步位移的QR 方法计算矩阵A (n-1)的全部特征值，也是A 的全部特征值，具体算法如下： 1. 给定精度水平0>ε和迭代最大次数L 。 2. 记n n ij n a A A ?-==][) 1()1()1(，令n m k ==,1。

数值分析大作业QR分解

题目:设计求取n n ?实数矩阵A 的所有特征值及其特征向量的数值算法，并以矩阵 20010-1-24A=0-2131 43 1?? ? ? ? ??? 为例进行具体的求解计算。 一、 算法分析： 一般而言，求取实数矩阵所有特征值的方法有雅克比法和QR 分解法，两者都是变换法。其中雅克比法只能求解对称矩阵的全部特征值和特征向量，而QR 则可用于更一般的矩阵特征值的求解，结合反幂法可进而求出各个特征向量。 二、 算法设计： 1、 原始实矩阵A R n n ?∈拟上三角化 为了减少求特征值和特征向量过程中的计算量，对生成的矩阵A 进行拟上三角化，得到拟上三角化矩阵A ’ 记A (1)=A ，并记A (r)的第r 列到第n 列的元素(1,2,...,;,1,...,)r ij a i n j r r n ==+。 对于r=1,2,…,n -2执行 (1) 若() (2,3,...,)r ir a i r r n =++全为零，则令A (r+1)= A (r)，转(5)；否则转(2)。 (2) 计算 令 ()2 ()() 1,1,1,sgn(0,sgn()=1) r r r r r r r r r r r c a a a ρ+++=-=，（若则取 (3) 令-0=r n r u u ?? ? ?? ，()()()-1,2,1(,,...,)r r r T n r r r r r r nr u a c a a ρ++=- (4) 计算 (r)(r)(r)T n-r r+1,r r r+2,r nr r T n-r T n-r n-r n-r n-r r+1r 1u = (a -c ,a ,...,a )ρ I H =I -2uu =H H =I -2u u A =HA H ?? ? ?? (5) 继续 算法执行完后，就得到与原矩阵A 相似的拟上三角矩阵A (n-1)。 2、 拟上三角矩阵QR 分解的求原矩阵的全部特征值 记A k 是对拟上三角矩阵进行QR 算法，产生的矩阵序列，A 0是起始拟上三角矩阵，

数值分析大作业 三、四、五、六、七

大作业 三 1. 给定初值0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件：fu.m function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件：dfu.m function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序Newton.m clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

while flag1==1 && m<=10^3 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)=ep flag=0; end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =- = ?-???解：Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while nerrorlim n=n+1; else break ; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数 非线性函数f function y=f(x) y=log((513+0.6651*x)/(513-0.6651*x))-x/(1400*0.0918); end

哈尔滨工程大学数值分析大作业2014-附fortran程序

B班大作业要求： 1. 使用统一封皮； 2. 上交大作业内容包含： 一摘要 二数学原理 三程序设计（必须对输入变量、输出变量进行说明；编程无语言要求，但程序要求通过）四结果分析和讨论 五完成题目的体会与收获 3. 提交大作业的时间：本学期最后一次课，或考前答疑；过期不计入成绩； 4. 提交方式：打印版一份；或手写大作业，但必须使用A4纸。 5. 撰写的程序需打印出来作为附录。

课程设计 课程名称： 设计题目： 学号： 姓名： 完成时间：

题目一：非线性方程求根 一 摘要 非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。本实验通过使用常用的求解方法二分法和Newton 法及改进的Newton 法处理几个题目，分析并总结不同方法处理问题的优缺点。观察迭代次数，收敛速度及初值选取对迭代的影响。 用Newton 法计算下列方程 （1） 310x x --= ， 初值分别为01x =，00.45x =，00.65x =； （2） 32943892940x x x +-+= 其三个根分别为1,3,98-。当选择初值02x =时给出结果并分析现 象,当6 510ε-=?，迭代停止。 二 数学原理 对于方程f(x)=0，如果f(x)是线性函数，则它的求根是很容易的。牛顿迭代法实质上是一种线性化方法，其基本思想是将非线性方程f(x)=0逐步归结为某种线性方程来求解。 设已知方程f(x)=0有近似根x k (假定k f'(x )0≠) ，将函数f(x)在点x k 进行泰勒展开，有 k k k f(x)f(x )+f'(x )(x-x )+≈??? 于是方程f(x)=0可近似的表示为 k k k f(x )+f'(x )(x-x )=0 这是个线性方程，记其根为x k+1，则x k+1的计算公式为 k+1k () x =x -'() k k f x f x ，k=0,1,2,… 这就是牛顿迭代法或简称牛顿法。

北航数值分析报告大作业第八题

北京航空航天大学 数值分析大作业八 学院名称自动化 专业方向控制工程 学号 学生姓名许阳 教师孙玉泉 日期2014 年11月26 日

一．题目 关于x , y , t , u , v , w 的方程组(A.3) ???? ?? ?=-+++=-+++=-+++=-+++79 .0sin 5.074.3cos 5.007.1cos sin 5.067.2cos 5.0y w v u t x w v u t y w v u t x w v u t (A.3) 以及关于z , t , u 的二维数表（见表A-1）确定了一个二元函数z =f (x , y )。 表A-1 二维数表 t z u 0 0.4 0.8 1.2 1.6 2 0 -0.5 -0.34 0.14 0.94 2.06 3.5 0.2 -0.42 -0.5 -0.26 0.3 1.18 2.38 0.4 -0.18 -0.5 -0.5 -0.18 0.46 1.42 0.6 0.22 -0.34 -0.58 -0.5 -0.1 0.62 0.8 0.78 -0.02 -0.5 -0.66 -0.5 -0.02 1.0 1.5 0.46 -0.26 -0.66 -0.74 -0.5 1. 试用数值方法求出f (x , y ) 在区域}5.15.0,8.00|), {≤≤≤≤=y x y x D （上的近似表达式 ∑∑===k i k j s r rs y x c y x p 00 ),( 要求p (x , y )以最小的k 值达到以下的精度 ∑∑==-≤-=10020 7210)],(),([i j i i i i y x p y x f σ 其中j y i x i i 05.05.0,08.0+==。 2. 计算),(),,(* ***j i j i y x p y x f (i =1,2,…,8 ; j =1,2,…,5) 的值，以观察p (x , y ) 逼 近f (x , y )的效果，其中j y i x j i 2.05.0,1.0**+==。

北航数值分析第二次大作业--QR分解

《数值分析A》

一、算法设计方案 整个程序主要分为四个函数，主函数，拟上三角化函数，QR分解函数以及使用双步位移求解矩阵特征值、特征向量的函数。因为在最后一个函数中也存在QR分解，所以我没有采用参考书上把矩阵M进行的QR分解与矩阵Ak的迭代合并的方法，而是在该函数中调用了QR分解函数，这样增强了代码的复用性，减少了程序长度；但由于时间关系，对阵中方法的运算速度没有进行深入研究。 1.为了减少QR分解法应用时的迭代次数，首先对给定矩阵进行拟上三角化处理。 2.对经过拟上三角化处理的矩阵进行QR分解。 3.注意到计算特征值与特征向量的过程首先要应用前面两个函数，于是在拟上三角化矩阵的基础上对QR分解函数进行了调用。计算过程中，没有采用goto语句，而是根据流程图采用其他循环方式完成了设计，通过对迭代过程的合并，简化了程序的循环次数，最后在计算特征向量的时候采用了列主元高斯消去法。

二、源程序代码 #include #include #include int i,j,k,l,m; //定义外部变量double d,h,b,c,t,s; double A[10][10],AA[10][10],R[10][10],Q[10][10],RQ[10][10]; double X[10][10],Y[10][10],Qt[10][10],M[10][10]; double U[10],P[10],T[10],W[10],Re[10]={0},Im[10]={0}; double epsilon=1e-12; void main() { void Quasiuppertriangular(double A[][10]); void QRdecomposition(double A[][10]); void DoublestepsQR(double A[][10]); int i,j; for(i=0;i<10;i++) { for(j=0;j<10;j++) { A[i][j]=sin(0.5*(i+1)+0.2*(j+1)); Q[i][j]=0; AA[i][j]=A[i][j]; } A[i][i]=1.5*cos(2.2*(i+1)); AA[i][i]=A[i][i];

数值计算大作业

数值计算大作业 题目一、非线性方程求根 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 +的计算公式为

数值分析作业答案.doc

第2章 插值法 1、当x=1，-1，2时，f(x)=0，-3，4，求f(x)的二次插值多项式。 （1）用单项式基底。 （2）用Lagrange 插值基底。 （3）用Newton 基底。 证明三种方法得到的多项式是相同的。 解：（1）用单项式基底 设多项式为:2 210)(x a x a a x P ++=， 所以：64 211111 1111122 2 211 200 -=-==x x x x x x A 3 76144 211111114241 13110111)() ()(22 221120 022 2 22 11 120 00-=-= ---==x x x x x x x x x f x x x f x x x f a 2 3694211111114411 31101111)(1)(1 )(122 221120 02 2 22112 001=--= --==x x x x x x x x f x x f x x f a 6 5654 2 1 1111114 2 1 3 11011111) (1)(1)(122 2 21120 022 11 00 2=--= ---==x x x x x x x f x x f x x f x a 所以f(x)的二次插值多项式为：26 52337)(x x x P ++-= （2）用Lagrange 插值基底 )21)(11() 2)(1())(())(()(2010210-+-+=----=x x x x x x x x x x x l )21)(11() 2)(1())(())(()(2101201------=----=x x x x x x x x x x x l ) 12)(12() 1)(1())(())(()(1202102+-+-=----= x x x x x x x x x x x l

数值分析大作业三四五六七完整版

数值分析大作业三四五 六七 HEN system office room 【HEN16H-HENS2AHENS8Q8-HENH1688】

大 作 业 三 1. 给定初值 0x 及容许误差 ，编制牛顿法解方程f (x )=0的通用 程序. 解：Matlab 程序如下： 函数m 文件： function Fu=fu(x) Fu=x^3/3-x; end 函数m 文件： function Fu=dfu(x) Fu=x^2-1; end 用Newton 法求根的通用程序 clear; x0=input('请输入初值x0:'); ep=input('请输入容许误差:'); flag=1; while flag==1 x1=x0-fu(x0)/dfu(x0); if abs(x1-x0)

end m=m+1; x0=x1; end if flag1==1||abs(x0)>=ep flag=0; end end fprintf('最大的sigma 值为：%f\n',sigma); 2．求下列方程的非零根 5130.6651()ln 05130.665114000.0918 x x f x x +?? =- = ?-???解：Matlab 程序为： （1）主程序 clear clc format long x0=765; N=100; errorlim=10^(-5); x=x0-f(x0)/subs(df(),x0); n=1; while nerrorlim n=n+1; else break ; end x0=x; end disp(['迭代次数: n=',num2str(n)]) disp(['所求非零根: 正根x1=',num2str(x),' 负根x2=',num2str(-x)]) （2）子函数 非线性函数f function y=f(x) y=log((513+*x)/*x))-x/(1400*; end （3）子函数 非线性函数的一阶导数df function y=df() syms x1 y=log((513+*x1)/*x1))-x1/(1400*; y=diff(y);