1、计算机在材料科学中的应用上机指导书计算机在材料科学中的应用上机指导书山东科技大学材料学院目录实验一 数值计算程序设计 4学时 3实验二 最小二乘法曲线拟合及回归分析 4学时 13实验三 相图热力学计算 4学时 16实验四 互联网在材料科学中的应用 2学时 18实验一 数值计算程序设计 4学时一、实验目的及任务要求1、实验目的：1）熟悉程序设计语言，了解简单的程序设计方法；2）熟悉迭代法求方程根及牛顿法求方程的算法设计；4）熟悉雅可比迭代法及其算法设计；5）熟悉追赶法解线性方程组及其算法设计。2、任务要求：1）熟悉机房环境；2）理解算法原理及程序编译过程。二、实验操作过程（一）方程求根 1学时1
2、、熟悉C语言编译软件界面回顾C语言程序的编写方法，熟悉C编译软件的特点。2、熟悉算法1）二分法求方程根基本思路： 假设方程 f(x)=0在区间a，b上有且仅有一根x*。（1）首先取区间的中点x0，则x0= (a+b) a,b对分为两个区间a,x0与x0,b；（2）判断：f(a)*f(x0)0：若不等式成立，则方程的根在x0,b，若不等式不成立，则方程的根在a,x0，根所在的新区间的端点仍然命名为a,b；（3）判定根是否达到预定的精度要求|xk+1-xk|，xk+1即可作为方程的根。2）牛顿法求根牛顿迭代法求根的方法是，将已知方程式f(x)=0看作函数式yf(x)。对函数f(x)在xk处进行一阶
3、泰勒展开,即可得:f(x)=f(xk)+f(xk)(x-xk)+函数f(x)=0时x值即为方程的根，即函数曲线与x轴的交点值。f(x)=f(xk)+f(xk)(x-xk)0x=xk-f(xk)/f(xk)设x=xk+1,则可得牛顿迭代求根的计算公式为：xk+1=xk-xk-f(xk)/f(xk)例如：已知方程式 f(x)=x3+4x2-10=0求得导函数 f(x)=3x2+8x因此可得牛顿迭代式为：xk+1=xk-几何解释：x1=x0-牛顿迭代法求根计算步骤：1）输入已知条件x0、N：x0初始迭代点，可判断其迭代收敛性；计算精度，一般取=0.00001N迭代预定最大次数，即当迭代次数超过N次之
4、后，则主观认为迭代发散或失败，此时停止计算。一般取N=10002）计算函数值f(x0)与导数值f(x0)3）迭代求解：x1=x0-4)精度判断：x1-x0|=，条件成立时，则输出计算结果；条件不成立，则x0=x1,k=k+1，返回第2步重新计算。3、程序设计1）读下面的程序，并改写采用二分法求解方程x3-x-1=0在区间0，2中的根（精度为e=10-6）二分法求一元n次方程的近似根:#include #include main() float x0,x1,x2,fx0,fx1,fx2; do printf(Enter x1,x2:); scanf(%f,%f,&x1,&x2); fx1=2*x
5、1*x1*x1-4*x1*x1+3*x1-6; fx2=2*x2*x2*x2-4*x2*x2+3*x2-6; while(fx1*fx20); /* 判断x1,x2 之间是否有根*/ do x0=(x1+x2)/2; /*取中点*/ fx0=2*x0*x0*x0-4*x0*x0+3*x0-6; if (fx0*fx1=1e-5); printf(%6.2f,x0); getch();2）读懂下面的牛顿法求方程的跟的程序，并改写求x3+4x2-10=0的根。给定条件：x0=1.0, =0.00001, N=1000计算结果：x*=1.3652300#includefloat solution(f
6、loat x)float x1,y,k;dok=6*x*x-8*x+3;y=2*x*x*x-4*x*x+3*x-6;x1=x-y/k;x=x1;while(fabs(y)0.001);return x;void main()float x;x=1.5;x=solution(x);printf(%fn,x);getch();（二）线性方程组求解 2学时1、熟悉算法1）雅可比迭代法算法设计线性方程组的矩阵表达式为AX=b，其一般形式为：Aijxj=bi (i=1,2, n j=1,2, n)将此式改写为迭代格式： （i=1, 2, ., n）迭代展开式为： x1=(b1-0-a12x2-a13x3
7、-a1nxn)/a11x2=(b2-a21x1-0-a23x3-a2nxn)/a22x3=(b3-a31x1-a32x2-0-a1nxn)/a33.xn=(bn-an1x1-an2x2-an3x3-0)/ann当给定一组x的初始迭代值时：， 按给定迭代式反复进行迭代求解即可得：， ， ， 如果迭代求解过程是收敛的， （i=1, 2, , n）雅可比迭代计算公式：（i=1, 2, , n）, (k=0, 1, 2, , )迭代求解的终止条件，也就是迭代求解的精度判断条件，即是用相邻两次解之差来决定，其判断条件有两种形式：1)取相邻两次解的每一个分量之差的绝对值|xi|，即：| （i=1, 2,
8、3, , n）2)取相邻两次解的每一个分量之差的绝对值中最大的一个|xi|max，即：|max （i=1, 2, 3, , n）2）追赶法求解三对角线方程组：三对角线方程组的形式如下：（1）追的过程即是消元过程： “a变0，b变1， c、d变成u和y即表示为：计算公式推导过程如下：第一式除b1得 x1+u1x2=y1其中 u1=c1/b1 y1=d1/b1第二式 第二式a2第一式，得：（b2-a2u1）x2+c2x3=d2-y1a2 x2+u2x3=y2其中：u2=c2/(b2-a2u1) y2=(d2-y1a2)/(b2-a2u1)第三式 第三式a3第二式，得：(b3-a3u2)x3+c3x
9、4=d3-a3y2 x3+u3x4=y3其中 u3=c3/(b3-a3u2) y3=(d3-a3y2)/(b3-a3u2)第n式 第n式an第n-1式，得：(bn-anun-1)xn=(dn-anyn-1) xn=yn其中 yn=(dn-anyn-1)/(bn-anun-1)经整理可得到下列方程组：u1=c1/b1, y1=d1/b1, ui=ci/(bi-aiui-1) (i=2,3, , n-1) yi=(di-aiyi-1)/(bi-aiui-1) (i=2, 3, .n)通过计算u1,u2,.,un-1及y1,y2,.yn，即可得到每个方程只含两个未知量的式子，而最后一式只含一个未知量
10、xn。因此，采用回代的方法即可求得全部解。由于计算出ui之后，ci就不再用了，因此，ui就用ci代替。同理，计算出yi之后，di就不再用了，因此，yi就用di代替：c1=c1/b1 d1=d1/b1t=bi-aici-1 ci=ci/t di=(di-aidi-1)/t (i=2, 3, , n-1)dn=(dn-andn-1)/(bn-ancn-1) （2）回代首先求 xn=yn在回代求 xn-1=yn-1-un-1xnxn-2=yn-2-un-2xn-1.x1=y1-u1x2因此，可得回代通式：xn=ynxi=yi-uixi+1 (i=n-1, n-2, , 1)由于在回代过程中，求得xi
11、之后，原来的yi不再用了，因此，xi可用yi代替。在消元过程中yi用di代替，ui用ci代替，所以xi也就用di代替了。这样，在全部计算公式中，ci代替了ui，di代替了yi、xi，因此，程序中仅只见到ci与di。回代公式可改写为：dn=dn di=di-cidi+1 (i=n-1, n-2, ., 1)。作业2: 计算例题：2、程序设计：1）读下面的方程，改写程序并以雅可比法求下列方程组：给定初始值 x(0)=0,0,0T选代求解为 x*=3,2,1T方程组： 5x1+2x2+x3=8 2x1+8x2-3x3=21 x1-3x2-6x3=1#include#include #include
12、#include #define EPS 1e-6/*允许误差*/ #define MAX 100/*迭代次数的最大值*/ float *Jacobi(float a34,int n); main() int i; float a34=5,2,1,8,2,8,-3,21,1,-3,-6,1;/*方程的系数*/ float *x; x=(float *)malloc(3*sizeof(float);/*动态申请内存，用于保存方程的解*/ x=Jacobi(a,3);/*调用雅可比函数*/ for(i=0;i3;i+)/*输出方程的解*/ printf(x%d=%ft,i,xi); printf(
13、n); getch(); float *Jacobi(float a34,int n) float *x,*y; float epsilon,s; int i,j,k=0; x=(float *)malloc(n*sizeof(float); y=(float *)malloc(n*sizeof(float); for(i=0;in;i+) xi=0; while(1) epsilon=0;/*容允误差*/ k+;/*迭代次数计数*/ for(i=0;in;i+) s=0; for(j=0;jn;j+) if(j=i) continue; s+=aij*xj; yi=(ain-s)/aii;/
14、*计算x(k+1)*/ epsilon+=(float)fabs(yi-xi); if(epsilon=MAX) printf(方程不收敛n); return y; for(i=0;in;i+) xi=yi; 2）读下面程序并分析： /追赶法解三对角阵方程组Ax=b /n-方程组的阶数。 /a1,a2,a3依次为三对角阵的三对角。 /成功返回1，失败返回0。#include main() int i,t;int RunEq(int n,double *a1,double *a2,double *a3,double *b,double *x); double a15=0,-1,-1,-1,-1;
15、double a25=4,4,4,4,4;double a35=-1,-1,-1,-1;double b5=100,200,200,200,100;double x5; t= RunEq(5,a1,a2,a3,b,x);if(t=1) for(i=0;i5;i+) printf(x%d=%fn,i,xi); getch(); int RunEq(int n,double *a1,double *a2,double *a3,double *b,double *x) int i; double l; for(i=1;in;i+) if(fabs(a2i-1)0.00000000001) retur
16、n 0; l=a1i/a2i-1; a2i-=a3i-1*l; bi-=bi-1*l; if(fabs(a2n-1)=0;i-) xi=(bi-a3i*xi+1)/a2i; return 1; 注意：a1的第一个元素需要空出来，因为矩阵第一行不含第一各对角的元素，同样a3的最后一个元素也要空出来。（二）excel在求解线性方程组的应用 1学时1、熟悉excel操作，并学会用excel求线性方程组的解1）方法1 x1+x2+2x3+3x4=13x1-x2-x3-2x4=-4 2x1+3x2-x3-x4=-6 x1+2x2+3x3-x4=-41.打开Excel。 2.由于在本方程组中未知数有4个，
17、所以预留4个可变单元格的位置A1-A4。 3.将活动单元格移至B1处，从键盘键入：=A1A22*A33*A4； 然后回车。 4.在B2处从键盘键入：=3 * A1A2A32 * A4；然后回车。 5.在B3处从键盘键入：=2 * A13 * A2A3A4；然后回车。 6.在B4处从键盘键入：=A12*A23*A3A4。 7.点击工具-规划求解，出现规划求解参数对话框。 8.对话框中第一栏为：设置目标单元格，在相应的框中填入$B$1。 9.对话框中第二栏为：等于；后有三个选项，依次为最大值，最小值，值为。根据题意B1表示方程组中第一个方程等号左边的表达式，它的值应为1，因此点击值为前的圆圈，输入
18、1。 10.对话框中第三栏为：可变单元格；我们预留的可变单元格为A1-A4，所以在可变单元格框内键入 $A $1：$A$4。 11.对话框中最后一栏为：约束；首先点击添加按钮，屏幕出现添加约束对话框。 12.在添加约束对话框的单元格引用位置键入：$B$2；在中间的下拉式菜单中选取=；在约束值处键入：4；然后按添加按钮，屏幕出现空白的添加约束对话框。 13.在添加约束对话框的单元格引用位置键入：$B$3；在中间的下拉式菜单中选取=；在约束值处键入：6；然后按添加按钮，出现空白的添加约束对话框。14.在添加约束对话框的单元格引用位置键入：$B$4；在中间的下拉式菜单中选取=；在约束值处键入：4；然
19、后按确定键，返回规划求解参数对话框。15.按求解键，出现求解结果对话框。此时在A1A4的位置依次为：1，1，0，1；这就是说，原方程组的解为：X1=1,X2=1,X3=0,X4=1。这样我们就求出了方程组的解。2）方法2理论基础：N阶线性方程组： AX = B 其中，A =aij ，i，j=1，2，n；B=(b1，b2， bn) ，为n维常数列向量。求：X=X1， X2，Xn；若系数行列式| A |0，则方程组必有唯一的解：X=A-1B。2X1+3X2+2X3+3X4=03X1+2X2- 2X3+3X4= 33X1+3X2+3X3- 4X4= 142X1- 2X2- 3X3- 3X4=7 1)
20、在Excel中输入系数方阵，如输入A1:D4。 2)判断线性方程组是否有解。 选择另外一个 单元元格，如E1，单击“常用” 具栏中“fx函数”按钮在“函数分类”中选择“数学与三角 数”类，然后选择“MDETERM”函数。在“Array”输入框中输入域A1:D4 。 3)求系数矩阵的逆A-1再选4行4列的一个区域F1:I4，单击“常用”工具栏中“fx函数”按钮在“函数分类” 中选择“数学与三角函数”类，然后选择“MINVERSE”函数。在“Array”输入框中输入区域A1：D4 并单击“确定”。将光标定位在编辑栏中所输入公式的结尾处，然后同时按下Ctrl，Shift，Enter 3个键，则在区域
21、F1:I4中显示出矩阵A 的逆矩阵A-1 的系数。4)求线性方程组的解再选定一个4行1列的区域如J1：J4 ，将列向量B输入到该区域中去。另外选择一个4行1列的区域如K1：K4，单击“常用”工具栏中“fx函数”按钮，选择“MMULT”函数。在“Array1 ”输入框中输入矩阵的逆A 所在区域F1 ：I4；在“Array2”输入框中输入列向量 所在的区域J1：J4，然后单击“确定”。将光标定位在编辑栏中所输入公式的结尾处，然后按下Ctrl，Shift，Enter 3个键，则区域K1：K4中显示出两个矩阵乘积结果，即方程的解。实验二 最小二乘法曲线拟合及回归分析 4学时一、实验目的及任务要求1、实
22、验目的：1）熟悉最小二乘法曲线拟合算法设计；2）熟悉该算法程序编写及调试方法。2、任务要求：1）熟悉一元及多元线性回归数值解法；2）熟悉用excel及origin进行回归分析的基本方法。二、实验操作过程1、一元线性回归（1）某种铝合金的含铝量为x(),其熔解温度为y（），由实验测得x与y地数据如下表所示。试用最小二乘法建立x与y的经验公式。2、多元线性回归P17 例2分析：求y对x1, x2, x3的线性回归方程。表22 去碳量y与天然矿石和烧结矿石的加入量x1, x2 及熔化时间x3实测数据 n=49，结合此题，写出多元线性拟合的基本步骤。3、求一个经验函数y=aebx ，a，b为常数，拟合
23、以下数据：x12345678y15.320.527.436.649.165.687.8117.6采用excel进行曲线拟合。4、 根据长期试验总结和金属材料理论，提出了耐热钢断裂时间与温度、持久强度的回归模型如下：其中y为断裂时间(h)；x为持久强度；T为试验温度(K)，R为气体常数。为求出25Cr2Mo1V钢的该回归模型，进行了27次试验，见表1。求该模型在给出条件下的系数，并对在工作温度为550和设计寿命为10万小时的条件下，此种耐热钢的持久强度做出估计。表1 25Cr2Mo1V耐热钢持久强度试验数据 No 温度/K 持久强度/MPa 断裂时间/h No 试验温度/K持久强度/ MPa 断
24、裂时间/h 1 823 400 113.5 15 853 270 937 2 823 380 163.5 16 853 250 1206.7 3 823 370 340.6 17 853 200 2044.6 4 823 360 561 18 873 300 182.2 5 823 350 953.8 19 873 270 350.7 6 823 350 1263.8 20 873 250 489.0 7 823 330 1902.8 21 873 200 958.7 8 823 310 2271.3 22 893 270 79.4 9 823 310 2466.5 23 893 250 15
25、0.4 10 823 270 3674.8 24 893 200 411.0 11 823 250 6368.7 25 893 150 1001.8 12 823 200 13862.0 26 893 120 1544.8 13 853 350 207.7 27 893 110 1795.0 14 853300 621.9结合课件，用origin进行回归分析。实验三 相图热力学计算 4学时一、实验目的及任务要求1、实验目的：1）熟悉在线相图查询网站；2）熟悉相图的热力学计算方法；3）熟悉简单相图计算程序。4）熟练掌握用Jade+origin软件分析XRD实验结果的方法2、任务要求：1）熟悉相图
26、热力学数值计算方法。2）熟悉Jade+origin的使用方法。二、实验操作过程1、在线相图计算： 2h1）FACT相图查询FACT是世界排名数一数二的热力学计算软件，提供在线热力学据计算，二元三元相图计算等等强大功能。下面的网址提供在线常用相图查询 http:/www.crct.polymtl.ca/FACT/documentation/试访问该网站并查询相图，在线查询Al2O3SiO2、Ni-Al、Fe-Al相图并下载。2）进入Knovel 电子图书数据库，进行相图查询和文献阅读。2、NiOMgO完全固溶体相图绘制NiOMgO为液、固相连续互溶二元体系，液相和固相均为理想溶液。已知NiO和MgO的熔点分别为1960和2800，熔化热分别为52.3kJ/Mol和77.4kJ/Mol。以纯液态NiO作为NiO的标准态，纯固态MgO作为MgO的标准态，则和近似计算公式为：=77400(1-=52300(1-设1为NiO，2为MgO，则有又因为11所以：1=(1) 求得：即可计算出NiOMgO完全固溶相图。读C源程序并写出各个模块的功能，尝试着改动程序，并运行看看对实验结果的影响。3、Jade标定物相方法-
