材料計(jì)算及模擬2

1、第二章 材料計(jì)算示例 n2.1 方程的求解與簡(jiǎn)單計(jì)算問(wèn)題 n2.2 線性回歸問(wèn)題 例2.1.1 求方程的解y=x3+10*x2-2*sin(x)-50分析：用matlab求解方程的符號(hào)解，所用函數(shù)為SOLVE(eqn1,eqn2,.,eqnN)SOLVE(eqn1,eqn2,.,eqnN,var1,var2,.,varN)SOLVE(eqn1,eqn2,.,eqnN,var1,var2,.varN)求解非線性方程的數(shù)值解(最小二乘法)所用函數(shù)fsolve函數(shù)文件格式(myfun為函數(shù)文件名) x = fsolve(myfun,x0,option)內(nèi)聯(lián)函數(shù)格式(fun為內(nèi)聯(lián)函數(shù)名) x = fs

2、olve(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,-10)或：fsolve(

3、ff,1)fsolve(ff,-1)fsolve(ff,-10)其中ff為函數(shù)文件，格式為：function f=ff(x)f=x3+10*x2-2*sin(x)-50;例例2.1.2 熱缺陷濃度的計(jì)算 熱缺陷是由于熱起伏引起的，并與溫度有關(guān)。故在某一溫度下，熱缺陷的數(shù)目可以用熱力學(xué)中自由能最小原理來(lái)進(jìn)行計(jì)算。熱缺陷濃度與溫度的關(guān)系式為 式中：N表示完整的單質(zhì)晶體的原子數(shù)目；n表示熱振動(dòng)形成的空位數(shù)；n/N表示熱缺陷在總結(jié)點(diǎn)中所占分?jǐn)?shù)，既熱缺陷濃度。Gf分別代表空位形成自由能或填充缺陷形成自由能。計(jì)算Gf分別等于1eV, 2eV, 4eV, 6eV和8eV時(shí)，溫度在100、1200和2000下

4、的缺陷濃度。Boltzmann常數(shù)k = 1.3810-23 J/K。1ev = 1.5910-19 J。 kTGNnf2exp例2.1.2 熱缺陷濃度的計(jì)算 例2.1.2 熱缺陷濃度的計(jì)算 例 2.1.3 自由能溫度關(guān)系式 CaCO3 分解反應(yīng)的自由能G與溫度的關(guān)系為:G=186.08+10.710-3TlnT+4.18710-6T2-5.23102T-1-0.245T(1) 計(jì)算溫度區(qū)間8001400K范圍內(nèi)的G。(間隔100K) (2) 繪出G與溫度的關(guān)系圖。 (3) 計(jì)算G=0的溫度條件。 例2.1.3 自由能溫度關(guān)系式 例2.1.4 微分方程的數(shù)值解 n常用函數(shù)有：ode23 精度1

5、0-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(t,y) tspan 變量

6、的間隔 y0 初始條件例2.1.4：解方程y = -y2+x, y(0)=1函數(shù)文件：(函數(shù)名：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程序定義函數(shù)：function dy = fode2(x,y)dy = zeros(2,1); % a column vectordy(1)

7、= 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 線性回歸問(wèn)題 例例2.2.1 實(shí)驗(yàn)測(cè)量含有不同量-氯代萘的一組線性聚乙烯試樣的熔點(diǎn)，得到的數(shù)據(jù)如下：- 氯代萘的體積分?jǐn)?shù)1 0.00 0.06 0.16 0.32 0.52 0.75 0.95聚乙烯的熔點(diǎn)Tm(K) 410.65 407.65 404.15 393.15 393.15 388.15 383.15假定熔點(diǎn)Tm與低分子稀釋劑的體積分?jǐn)?shù)1的關(guān)系為式中，u和1分別

8、是高分子和低分子稀釋劑的密度。R氣體常數(shù)。1是高分子和稀釋劑的相互作用參數(shù)。Hu是每摩爾重復(fù)單元的熔融熱。求回歸系數(shù)。)(1121111uuommHRTT2.2 線性回歸問(wèn)題 )(1121111uuommHRTT21111111uuuuommHRHRTT22110 xbxbby2.2 線性回歸問(wèn)題 mmxbxbxbby.22110nnmnnmmmmyyyyxxxxxxxxxxxxbbbb32121332312222111211210;1111;yXBYXX)(XBTT12.2 線性回歸問(wèn)題 Tm = 410.65 407.65 404.15 393.15 393.15 388.15 383.1

9、5;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 線性回歸問(wèn)題 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 r = corrcoef

10、(Y,Ytheory) % 相關(guān)系數(shù)x = 0: 0.01: 1;plot(fai, Y, o, x, polyval(P,x), -)例例2.2.2 用DSC法測(cè)試PET的結(jié)晶動(dòng)力學(xué)數(shù)據(jù)，結(jié)果如下表，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方程求解結(jié)晶動(dòng)力學(xué)參數(shù). Avrami方程可寫成：對(duì)于線性模型：用最小二乘法中的矩陣法( )計(jì)算出Avrami方程中的參數(shù)Z和n。

11、tnZlnln)1ln(lnaxbby10YXXXBTT1)(2.2 線性回歸問(wèn)題 例例2.2.3 熱膨脹問(wèn)題熱膨脹問(wèn)題 研究米尺基準(zhǔn)器的線膨脹系數(shù)，得出在不同溫度時(shí)該基準(zhǔn)器的長(zhǎng)度修正值可用如下公式表示式中，x為0 oC時(shí)米尺基準(zhǔn)器的修正值(單位mm)；y和z為溫度系數(shù)；T為溫度。在不同溫度時(shí)米尺基準(zhǔn)器長(zhǎng)度的修正值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.09 299.84試用線性回

12、歸法求出x、y、z的值。 2zTyTxL2.2 線性回歸問(wèn)題例2.2.4 Avrami方程的加權(quán)最小二乘法擬合 n用差示掃描量熱法(DSC)研究聚對(duì)苯二甲酸乙二酯在232.4的等溫結(jié)晶過(guò)程，由結(jié)晶放熱峰原始曲線獲得結(jié)晶時(shí)間t與相應(yīng)的相對(duì)平衡結(jié)晶度的數(shù)據(jù)如下： 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試用加權(quán)最小二乘法求出Avrami方程中的參數(shù)(n和Z) Avrami方程:)exp(1nZta線性形式:tnZlnln)1ln(lna模型:bxby0

13、YXXXBTT1)(矩陣解:加權(quán)最小二乘法矩陣解:WYXWXXBTT1)(函數(shù)隨機(jī)誤差傳遞公司求函數(shù)y的方差:22222)1ln()1 (1aaayyy的權(quán)函數(shù):22)1ln()1(11)(aaayW權(quán)函數(shù)中的常數(shù)項(xiàng)可省略,故:2)1ln()1()(aaaW權(quán)函數(shù)W的矩陣形式:nwww0021W對(duì)角線上的矩陣元素可根據(jù)權(quán)函數(shù)計(jì)算得到擬合效果的判據(jù):令擬合殘差平方和為 niiixbbywQ1210)(m為待定參數(shù)的個(gè)數(shù)。一元線性回歸m=2.mnQ2單位權(quán)方差為 方差協(xié)方差矩陣 122)()(WXXBsT相關(guān)系數(shù)的平方： niiiyywQr122)(1是 yi的加權(quán)平均值：niiniiiwywy1

14、1yclearclct = 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 = 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