材料計算與模擬

1、第二章 材料計算示例 n2.1 方程的求解與簡單計算問題 n2.2 線性回歸問題 n2.3 非線性最小二乘法例2.1.1 求方程的解y=x3+10*x2-2*sin(x)-50分析：用matlab求解方程的符號解，所用函數(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為

2、內(nèi)聯(lián)函數(shù)名) 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為函數(shù)文件，格式為：function f=ff(x)f=x3+10*x2-2*sin(x)-50;例例2.1.2 熱缺陷濃度的計算 熱缺陷是由于熱起伏引起的，并與溫度有關。故在某一溫度下，熱缺陷的數(shù)目可以用熱力學中自由能最小原理來進行計算。熱缺陷濃度與溫度的關系式為 式中：N表示完整的單質(zhì)晶體的原子數(shù)目；n表示熱振動形成的空位數(shù)；n/N表示熱缺陷在總結點中所占分數(shù)，既熱缺陷濃度。Gf分別代表空位形成自由能或填充缺陷形成自由能。計算Gf分別等于1eV, 2eV, 4eV, 6eV和8eV時，溫度在1

4、00、1200和2000下的缺陷濃度。Boltzmann常數(shù)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) 計算溫度區(qū)間8001400K范圍內(nèi)的G。(間隔100K) (2) 繪出G與溫度的關系圖。 (3) 計算G=0的溫度條件。 例2.1.3 自由能溫度關系式 例2.1.4 微分方程的數(shù)值解 n常用

5、函數(shù)有：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函數(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

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 實驗測量含有不同量-氯代萘的一組線性聚乙烯試樣的熔點，得到的數(shù)據(jù)如下：- 氯代萘的體積分數(shù)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與低分子稀釋劑的體積分數(shù)

8、1的關系為式中，u和1分別是高分子和低分子稀釋劑的密度。R氣體常數(shù)。1是高分子和稀釋劑的相互作用參數(shù)。Hu是每摩爾重復單元的熔融熱。求回歸系數(shù)。)(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) % 相關系數(shù)x = 0: 0.01: 1;plot(fai, Y, o, x, polyval(P,x), -)例例2.2.2 用DSC法測試PET的結晶動力學數(shù)據(jù)，結果如下表，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方程求解結晶動力學參數(shù). Avrami方程可寫成：對于線性模型：用最小二乘法中的矩陣法( )計算出Avr

11、ami方程中的參數(shù)Z和n。tnZlnln)1ln(lnaxbby10YXXXBTT1)(2.2 線性回歸問題 例例2.2.3 熱膨脹問題熱膨脹問題 研究米尺基準器的線膨脹系數(shù)，得出在不同溫度時該基準器的長度修正值可用如下公式表示式中，x為0 oC時米尺基準器的修正值(單位mm)；y和z為溫度系數(shù)；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與相應的相對平衡結晶度的數(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試用加權最小二乘法求出Avrami方程中的參數(shù)(n和Z) Avrami方程:)exp(1nZta線性形式:tnZlnln)1l

13、n(lna模型:bxby0YXXXBTT1)(矩陣解:加權最小二乘法矩陣解:WYXWXXBTT1)(函數(shù)隨機誤差傳遞公司求函數(shù)y的方差:22222)1ln()1 (1aaayyy的權函數(shù):22)1ln()1(11)(aaayW權函數(shù)中的常數(shù)項可省略,故:2)1ln()1()(aaaW權函數(shù)W的矩陣形式:nwww0021W對角線上的矩陣元素可根據(jù)權函數(shù)計算得到擬合效果的判據(jù):令擬合殘差平方和為 niiixbbywQ1210)(m為待定參數(shù)的個數(shù)。一元線性回歸m=2.mnQ2單位權方差為 方差協(xié)方差矩陣 122)()(WXXBsT相關系數(shù)的平方： 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的結晶動力學數(shù)據(jù)，結果如下表，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方程求解結晶動力學參數(shù). Avrami方程可寫成：用非線性最小二乘法計算出Avrami方程中的參數(shù)Z和n。分析：Matlab的非線性擬合函數(shù)：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 假設有一組實驗數(shù)據(jù)假設有一組實驗數(shù)據(jù)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