材料计算与模拟

1、第二章 材料计算示例 n2.1 方程的求解与简单计算问题 n2.2 线性回归问题 n2.3 非线性最小二乘法例2.1.1 求方程的解y=x3+10*x2-2*sin(x)-50分析：用matlab求解方程的符号解，所用函数为SOLVE(eqn1,eqn2,.,eqnN)SOLVE(eqn1,eqn2,.,eqnN,var1,var2,.,varN)SOLVE(eqn1,eqn2,.,eqnN,var1,var2,.varN)求解非线性方程的数值解(最小二乘法)所用函数fsolve函数文件格式(myfun为函数文件名) x = fsolve(myfun,x0,option)内联函数格式(fun为

2、内联函数名) x = fsolve(fun,x0,option);fplot(x3+10*x2-2*sin(x)-50,-15 5),grid on-15-10-505-1200-1000-800-600-400-2000200400f=inline(x3+10*x2-2*sin(x)-50,x);fplot(f,-15 5); grid onfsolve(f,-10)fsolve(f,-1)fsolve(f,1)或：fsolve(x3+10*x2-2*sin(x)-50,1)fsolve(x3+10*x2-2*sin(x)-50,-1)fsolve(x3+10*x2-2*sin(x)-50,

3、-10)或：fsolve(ff,1)fsolve(ff,-1)fsolve(ff,-10)其中ff为函数文件，格式为：function f=ff(x)f=x3+10*x2-2*sin(x)-50;例例2.1.2 热缺陷浓度的计算 热缺陷是由于热起伏引起的，并与温度有关。故在某一温度下，热缺陷的数目可以用热力学中自由能最小原理来进行计算。热缺陷浓度与温度的关系式为 式中：N表示完整的单质晶体的原子数目；n表示热振动形成的空位数；n/N表示热缺陷在总结点中所占分数，既热缺陷浓度。Gf分别代表空位形成自由能或填充缺陷形成自由能。计算Gf分别等于1eV, 2eV, 4eV, 6eV和8eV时，温度在1

4、00、1200和2000下的缺陷浓度。Boltzmann常数k = 1.3810-23 J/K。1ev = 1.5910-19 J。 kTGNnf2exp例2.1.2 热缺陷浓度的计算 例2.1.2 热缺陷浓度的计算 例 2.1.3 自由能温度关系式 CaCO3 分解反应的自由能G与温度的关系为:G=186.08+10.710-3TlnT+4.18710-6T2-5.23102T-1-0.245T(1) 计算温度区间8001400K范围内的G。(间隔100K) (2) 绘出G与温度的关系图。 (3) 计算G=0的温度条件。 例2.1.3 自由能温度关系式 例2.1.4 微分方程的数值解 n常用

5、函数有：ode23 精度10-3 , ode45 精度 10-6 n格式： T,Y = solver(odefun,tspan,y0)T,Y = solver(odefun,tspan,y0,options)T,Y = solver(odefun,tspan,y0,options,p1,p2.)T,Y,TE,YE,IE = solver(odefun,tspan,y0,options)sol = solver(odefun,t0 tf,y0.)where solver is one of ode45, ode23, ode113odefun方程: 例如y=f(t,y)或M(t,y)y = f(

6、t,y) tspan 变量的间隔 y0 初始条件例2.1.4：解方程y = -y2+x, y(0)=1函数文件：(函数名：fode1)function dy = fode1(x,y) dy=-y2+x主程序：x,y = ode23(fode1,0,10,1);plot(x,y)xlabel(x);ylabel(y);例2.1.5解方程： 可化成： 0)0(2)0()1(10002yyyyyy0)0(2)0()1(1000211221221yyyyyyyy例2.1.5：Matlab程序定义函数：function dy = fode2(x,y)dy = zeros(2,1); % a column

7、 vectordy(1) = y(2);dy(2) = 1000*(1 - y(1)2)*y(2) - y(1);主程序:x,y = ode23(fode2,0,5,2,0);plot(x,y(:,1)xlabel(x);ylabel(y); 2.2 线性回归问题 例例2.2.1 实验测量含有不同量-氯代萘的一组线性聚乙烯试样的熔点，得到的数据如下：- 氯代萘的体积分数1 0.00 0.06 0.16 0.32 0.52 0.75 0.95聚乙烯的熔点Tm(K) 410.65 407.65 404.15 393.15 393.15 388.15 383.15假定熔点Tm与低分子稀释剂的体积分数

8、1的关系为式中，u和1分别是高分子和低分子稀释剂的密度。R气体常数。1是高分子和稀释剂的相互作用参数。Hu是每摩尔重复单元的熔融热。求回归系数。)(1121111uuommHRTT2.2 线性回归问题 )(1121111uuommHRTT21111111uuuuommHRHRTT22110 xbxbby2.2 线性回归问题 nnxbxbxbby.22110mnmnnmmmnyyyyxxxxxxxxxxxxbbbb32121332312222111211210;1111;yXBYXX)(XBTT12.2 线性回归问题 Tm = 410.65 407.65 404.15 393.15 393.15

9、 388.15 383.15;fai = 0.00 0.06 0.16 0.32 0.52 0.75 0.95;X = ones(7,3);X(:,2) = fai;X(:,3) = fai.2;Y = 1./Tm;format short e2.2 线性回归问题 format short eB1 = (X*X)-1*X*Y % method 1B2 = (X*X)(X*Y) % method 2P,S = polyfit(fai,Y,2) % method 3B3 = XY %method 4Ytheory = polyval(P,fai); % Y的理论值 Yval = Ytheory,Y

10、 r = corrcoef(Y,Ytheory) % 相关系数x = 0: 0.01: 1;plot(fai, Y, o, x, polyval(P,x), -)例例2.2.2 用DSC法测试PET的结晶动力学数据，结果如下表，t(min) 0.288 0.534 0.822 1.110 1.398 1.686 1.972 2.342 2.712 3.206a 0.015 0.101 0.348 0.647 0.841 0.925 0.958 0.978 0.988 0.995用Avrami方程求解结晶动力学参数. Avrami方程可写成：对于线性模型：用最小二乘法中的矩阵法( )计算出Avr

11、ami方程中的参数Z和n。tnZlnln)1ln(lnaxbby10YXXXBTT1)(2.2 线性回归问题 例例2.2.3 热膨胀问题热膨胀问题 研究米尺基准器的线膨胀系数，得出在不同温度时该基准器的长度修正值可用如下公式表示式中，x为0 oC时米尺基准器的修正值(单位mm)；y和z为温度系数；T为温度。在不同温度时米尺基准器长度的修正值L如下表所示：T(oC) 0.551 5.363 10.459 14.277 17.806 22.103 24.633 28.986 34.417L(mm) 5.70 47.61 91.49 124.25 154.87 192.64 214.57 252.0

12、9 299.84试用线性回归法求出x、y、z的值。 2zTyTxL2.2 线性回归问题例2.2.4 Avrami方程的加权最小二乘法拟合 n用差示扫描量热法(DSC)研究聚对苯二甲酸乙二酯在232.4的等温结晶过程，由结晶放热峰原始曲线获得结晶时间t与相应的相对平衡结晶度的数据如下： t (min)7.611.317.421.625.6 27.6 31.635.636.638.1(%) 3.41 11.534.754.972.7 80.091.097.398.299.3试用加权最小二乘法求出Avrami方程中的参数(n和Z) Avrami方程:)exp(1nZta线性形式:tnZlnln)1l

13、n(lna模型:bxby0YXXXBTT1)(矩阵解:加权最小二乘法矩阵解:WYXWXXBTT1)(函数随机误差传递公司求函数y的方差:22222)1ln()1 (1aaayyy的权函数:22)1ln()1(11)(aaayW权函数中的常数项可省略,故:2)1ln()1()(aaaW权函数W的矩阵形式:nwww0021W对角线上的矩阵元素可根据权函数计算得到拟合效果的判据:令拟合残差平方和为 niiixbbywQ1210)(m为待定参数的个数。一元线性回归m=2.mnQ2单位权方差为 方差协方差矩阵 122)()(WXXBsT相关系数的平方： niiiyywQr122)(1是 yi的加权平均值

14、：niiniiiwywy11yclearclct = 7.6 11.3 17.4 21.6 25.6 27.6 31.6 35.6 36.6 38.1;alfa = 3.41 11.5 34.7 54.9 72.7 80.0 91.0 97.3 98.2 99.3;X(1:10,1) = 1;X(1:10,2) = log(t);Y=log(-log(1-alfa./100);for i = 1:10 w(i) = (1-alfa(i)/100)*log(1-alfa(i)/100)2; W(i,i) = w(i);endB = (X*W*X)-1*X*W*Y;Z = exp(B(1)n =

15、B(2)Q = 0;my = 0;for i = 1:10 Q=Q+w(i)*(Y(i)-(B(1)+B(2)*log(t(i)2; my = my + Y(i);endmy = my/10;Qm = 0;for i = 1:10 Qm = Qm + w(i)*(Y(i)-my)2;endr = 1 - Q/Qm Z = 1.0128e-004n = 2.9188r = 0.9994 plot(t,w,o);hold onxx = 7:0.01:38.5;yy = spline(t,w,xx);plot(xx,yy,-)51015202530354000.020.040.060.080.10.

16、120.14timeweighted function w(a)例例2.3.1 用DSC法测试PET的结晶动力学数据，结果如下表，t(min) 0.288 0.534 0.822 1.110 1.398 1.686 1.972 2.342 2.712 3.206a 0.015 0.101 0.348 0.647 0.841 0.925 0.958 0.978 0.988 0.995用Avrami方程求解结晶动力学参数. Avrami方程可写成：用非线性最小二乘法计算出Avrami方程中的参数Z和n。分析：Matlab的非线性拟合函数：lsqcurvefit(fun,x0,xdata,ydata

17、)exp(1nZta2.3 非线性曲线拟合 程序t=0.288 0.534 0.822 1.110 1.398 1.686 1.972 2.342 2.712 3.206;alfa=0.015 0.101 0.348 0.647 0.841 0.925 0.958 0.978 0.988 0.995;f=inline(1-exp(-a(1)*x.a(2),a,x)a=lsqcurvefit(f,1 1,t,alfa);a运行结果：a = 0.7349 2.7248例例2.3.2 假设有一组实验数据假设有一组实验数据x 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1y 2.3201 2.6470 2.9707 3.2885 3.6008 3.9090