中南大學材料專業(yè)MATLAB實踐4號

1、一、MATLAB程序設計實踐Matlab基礎班級：學號：姓名：表示多晶體材料織構的三維取向分布函數(shù)（ff（1，2）是一個非常復雜的函數(shù)，難以精確的用解析函數(shù)表達，通常采用離散空間函數(shù)值來表示取向分布函數(shù)，Data.txt是三維取向分布函數(shù)的一個實例。由于數(shù)據(jù)量非常大，不便于分析，需要借助圖形來分析。請你編寫一個matlab程序畫出如下的幾種圖形來分析其取向分布特征：（1）用Slice函數(shù)給出其整體分布特征；（2）用pcolor或contour函數(shù)分別給出(20, 5, 10, 15, 20, 25, 30, 35 90)切面上f分布情況(需要用到subplot函數(shù))；(3) 用plot函數(shù)給出

2、沿取向線(1=090，45，20)的f分布情況。備注：表示以下數(shù)據(jù)為20的數(shù)據(jù)，即f（1，0）此方向表示f隨著從0,5,10,15, 20 到90的變化而變化此方向表示f隨著1從0,5,10,15, 20 到90的變化而變化數(shù)據(jù)說明部分，與作圖無關data.txt數(shù)據(jù)格式說明將文件Data.txt內的數(shù)據(jù)按照要求讀取到矩陣f(phi1,phi,phi2)中，代碼如下：fid=fopen('data.txt'); %讀取數(shù)據(jù)文件Data.txtfor i=1:18 tline=fgetl(fid);endphi1=1;phi=1;phi2=1;line=0;f=zeros(19,

3、19,19);while feof(fid) tline=fgetl(fid); data=str2num(tline); line=line+1; if mod(line,20)=1 phi2=(data/5)+1; phi=1; elsefor phi1=1:19 f(phi1,phi,phi2)=data(phi1);endphi=phi+1;endendfclose(fid);建立新的腳本，將以上代碼保存為m文件，命名為readtext.m，在MATLAB運行后。其變量f結果為：(1). 用Slice函數(shù)給出其整體分布特征程序代碼：fopen('readtext.m')

4、;readtext;x,y,z=meshgrid(0:5:90,0:5:90,0:5:90); slice(x,y,z,f,45,90,45,90,0,45) %運用slice函數(shù)繪制圖形將程序代碼保存為code1_1.m文件，其運行結果：(2). 用subplot函數(shù)和pcolor函數(shù)給出(20, 5, 10, 15, 20, 25, 30, 35 90)切面上f分布情況；程序代碼：fopen('readtext.m'); readtext; for i=1:19 subplot(5,4,i) pcolor(f(:,:,i) %運用pcolor函數(shù)繪制圖形end將上述代碼保存

5、為code1_2_1.m文件，其運行結果：. 使用subplot函數(shù)和contour函數(shù)給出(20, 5, 10, 15, 20, 25, 30, 35 90)切面上f分布情況；程序代碼：fopen('readtext.m'); %運用contour函數(shù)繪制圖形 readtext; for i=1:19 subplot(5,4,i) contour(f(:,:,i) end將上述代碼保存為code1_2_2.m文件，運行結果：（3）. 用plot函數(shù)給出沿取向線(1=090，45，20)的f分布情況。程序代碼：fopen('readtext.m'); readt

6、ext; plot(0:5:90,f(:,10,1),'-bo') %運用plot函數(shù)繪制圖形 text(60,6,'phi=45 phi2=0')將上述代碼保存為code1_3.m文件，運行結果：二 MATLAB程序設計實踐科學計算（04）班級、編程實現(xiàn)以下科學計算算法，并舉一例應用之。（參考書籍精通科學計算，王正林等著，電子工業(yè)出版社，年）“切比雪夫逼近”算法說明：當一個連續(xù)函數(shù)定義在區(qū)間-1,1上時，它可以展開成切比雪夫級數(shù)。即：其中為次切比雪夫多項式，具體表達可通過遞推得出： ， 它們之間滿足如下的正交關系： 在實際應用中，可根據(jù)所需的精度來截取有限的項

7、數(shù)，切比雪夫級數(shù)中的系數(shù)由下式?jīng)Q定： 流程圖：結束開始定義sym型變量t以及切比雪夫多項式矩陣T(n)，并規(guī)定前兩項為1，t輸出k項切比雪夫多項式求出逼近x0處六位精度下的值根據(jù)前兩項多項式計算其系數(shù)c(1)，c(2)是輸入三項?否由第一項和第二項多項式及其系數(shù)計算一級切比雪夫逼近值是i>k+1？i從3漸漸增加到k+1 否 是 由T(i-1)和T(i-2)計算T(i) 由T(i)和c(i)計算i級切比雪夫逼近值。由T(i)計算其系數(shù)c(i)程序源代碼（m文件）：function f = Chebyshev(y,k,x0)%用切比雪夫多項式逼近已知函數(shù)%已知函數(shù)：y%逼近已知函數(shù)所需項數(shù)：

8、k%逼近點的x坐標：x0%求得的切比雪夫逼近多項式或在x0處的逼近值：f syms t;%定義sym型變量tT(1:k+1) =t;T(1) = sym('1');T(2) = t;%定義切比雪夫多項式矩陣T(n)，并規(guī)定前兩項為1，tc(1:k+1) = sym('0');c(1)=int(subs(y,findsym(sym(y),sym('t')*T(1)/sqrt(1-t2),t,-1,1)/pi;c(2)=2*int(subs(y,findsym(sym(y),sym('t')*T(2)/sqrt(1-t2),t,-1,

9、1)/pi;%根據(jù)前兩項多項式計算其系數(shù)c(1)，c(2)f = c(1)+c(2)*t;%由第一項和第二項多項式及其系數(shù)計算切比雪夫逼近值 for i=3:k+1 %從第3項到第k+1項進行循環(huán) T(i) = 2*t*T(i-1)-T(i-2);%計算第i項切比雪夫多項式矩陣T(i) c(i) = 2*int(subs(y,findsym(sym(y),sym('t')*T(i)/sqrt(1-t2),t,-1,1)/2; %計算第i項系數(shù) f = f + c(i)*T(i);%計算切比雪夫逼近值 f = vpa(f,6); if(i=k+1) if(nargin = 3)

10、f = subs(f,'t',x0);%輸入三項，則輸出x0處六級精度的切比雪夫逼近值 end f = vpa(f,6); endend根據(jù)代碼可知，在MATLAB中編程實現(xiàn)的切比雪夫逼近法函數(shù)為：Chebyshev。功能：用切比雪夫多項式逼近已知函數(shù)。調用格式：其中，y為已知函數(shù)； k為逼近已知函數(shù)所需項數(shù)； f是求得的切比雪夫逼近多項式或在x0處的逼近值。應用實例：切比雪夫應用實例。用切比雪夫公式（取6項）逼近函數(shù)，并求當x=1時的函數(shù)值。解：利用程序求解方程，在MATLAB命令窗口中輸入：>> Chebyshev('1/(2-x)',5) %調

11、用創(chuàng)建的函數(shù)，輸出切比雪夫多項式的5個項結果：ans =0.277013*t + 0.0647768*t*(2.0*t2 - 1.0) - 0.00501051*t*(2.0*t*(t - 2.0*t*(2.0*t2 - 1.0) + 2.0*t2 - 1.0) - 0.0186995*t*(t - 2.0*t*(2.0*t2 - 1.0) + 0.24175*t2 + 0.456475再在MATLAB命令窗口中輸入：>> Chebyshev('1/(2-x)',5,1) %調用創(chuàng)建的函數(shù)，輸出當x=1時的函數(shù)值逼近值結果：ans =1.06372在MATLAB命令

12、窗口計算符號表達式x=1的極限：>>limit(sym('1/(2-x)'),'x',1) %計算其真值，并與其進行比較結果：ans = 1輸出結果：絕對誤差：1.06372-1=0.06372；相對誤差：0.06372/1*100%=6.372%、編程解決以下科學計算問題。1)問題分析：1.相關參數(shù)字母設定：如圖，由在P的作用，點A移動到點A，設在水平方向移動距離為dx，在豎直方向移動距離為dy，P與水平的夾角為。設各桿的軸力為Ni(i=1,2,3,n)；各桿與水平面夾角分別i(i=1,2,3,n)；各桿的變化長度分別為dLi(i=1,2,3,n)

13、；2.解題思路：（1）取一根桿分析，由題意可知每根桿的伸長量為dx和dy在桿上的分量；故可得方程：dLi=NiLiEAi=dxcosi+dysini，i=1,2,3n (2)將各軸力和外力P水平（x）豎直（y）正交分解，由受力平衡可得：x方向有方程： i=1nNicosi=pcos y方向有方程： i=1nNisini=psin （3）聯(lián)立上述可得n+2等式，可解得n個軸力，以及A點位移dx和dy，聯(lián)立的線性方程組如下：.i=1nNicosi=Pcosi=1nNisini=Psin（4）建立Pcos,Psin，0，0，00的常數(shù)矩陣，以及如下所示的系數(shù)矩陣：N1cos1N2cos2N3cos3

14、Nncosn00N1sin1N2sin2N3sin3Nnsinn00N1L1EA100-cos1-sin10N2L2EA20-cos2-sin200N3L3EA30-cos3-sin3000NnLnEAn-cosn-sinn（5）再用求逆法求解此線性方程組，即用常數(shù)矩陣除以系數(shù)矩陣，得出結果。程序源代碼：Ei=input('請輸入各桿的剛度:(注意用括起來)'); %輸入剛度矩陣EiLi=input('請輸入各桿的長度:(注意用括起來) ');%輸入桿的長度矩陣LiAi=input('請輸入各桿的橫截面積:(注意用括起來) ');%輸入桿的橫截面

15、積矩陣Aiai=input('請輸入各桿與水平面的夾角:(注意用括起來) ');%輸入桿與水平面的夾角矩陣aiP=input('請輸入外力P: ');%輸入外加力Pa=input('請輸入P與水平面的夾角: ');%輸入外加力P與x的夾角n1=length(Ei);n2=length(Li);n3=length(Ai);n4=length(ai);if(n1=n2|n1=n3|n1=n4) disp('輸入數(shù)據(jù)錯誤')else n=n1; end%判斷賦值維度是否一致Ki=Li./(Ei.*Ai);C=zeros(n+2,1);C

16、(1,1)=P*cos(a);C(2,1)=P*sin(a);C(3:n+2,1)=zeros(n,1);%建立方程組等號右邊常數(shù)的矩陣D=zeros(n+2,n+2);D(1,:)=cos(ai),0,0;D(2,:)=sin(ai),0,0;for(i=1:n) D(i+2,i)=Ki(i);endD(3:n+2,n+1)=(-cos(ai);D(3:n+2,n+2)=(-sin(ai);%建立方程組系數(shù)矩陣x=DC;x=x'%求解該線性方程組，得出位移以及每根桿的軸力disp('節(jié)點在x、y方向上的位移分別：')x(n+1:n+2)disp('各桿的軸力分

17、別為：')x(1:n)%輸出結果流程圖：開始結束輸入所需數(shù)據(jù)：各桿的剛度Ei、橫截面積Ai、長度Li，以及夾角ai，外加力P以及其與x方向夾角輸出A點位移，以及各桿軸力輸出：輸入數(shù)據(jù)錯誤用求逆法解線性微分方程輸入Ei，Ai，Li，i維度是否一致否建立系數(shù)矩陣D，以及常數(shù)矩陣C是應用實例：設三根桿組成的支架如下圖所示，掛一重物P=3000N。設L=3m，各桿的橫截面積分別為：A1=150×10-6m2，A2=200×10-6m2，A3=300×10-6m2，材料的彈性模量均為E=200×109N/m2，求各桿所受力的大小以及C點位移解：將前面的程序

18、源代碼保存為main.m文件運行，根據(jù)提示依次輸入題中數(shù)據(jù)：200e9,200e9,200e9;3/sin(pi/3),3/sin(pi/2),3/sin(pi/4);150e-6,200e-6,300e-6;pi/3,pi/2,3*pi/4;3000;0;輸出結果：節(jié)點在x、y方向上的位移分別：ans = 1.0e-03 * 0.3399 0.0420各桿的軸力分別為：ans = 1.0e+03 * 1.7865 0.5595 -2.97942)開始流程圖： 輸入MATLAB指令，調用函數(shù)進行計算 程序出錯，返回修改 得到結果，結束 （1）0112exp（-x22）dx程序代碼：>> fun=inline('