一、《MATLAB程序设计实践》Matlab基础班级：学号：姓名：表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，Data.txt是三维取向分布函数的一个实例。

由于数据量非常大，不便于分析，需要借助图形来分析。

请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：（1）用Slice函数给出其整体分布特征；（2）用pcolor或contour函数分别给出(φ2＝0, 5, 10, 15, 20, 25,30, 35 …90)切面上f分布情况(需要用到subplot函数)；(3) 用plot函数给出沿α取向线(φ1=0~90，φ＝45，φ2＝0)的f分布情况。

流程图解：（1）将文件Data.txt内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下：fid=fopen('data.txt'); %读取数据文件Data.txtfor i=1:18tline=fgetl(fid);endphi1=1;phi=1;phi2=1;line=0;f=zeros(19,19,19);while ~feof(fid)tline=fgetl(fid);data=str2num(tline);line=line+1;if mod(line,20)==1phi2=(data/5)+1;phi=1;elsefor phi1=1:19f(phi1,phi,phi2)=data(phi1);endphi=phi+1;endendfclose(fid);将以上代码保存为readtext.m文件并在MATLAB中运行，运行结果如下图所示：fopen('readtext.m');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函数绘制图形运行结果如右图所示（2）将以下代码保存为code1_2_1.m文件：fopen('readtext.m');readtext;for i=1:19subplot(5,4,i)pcolor(f(:,:,i)) %运用pcolor函数绘制图形end运行结果如右图所示fopen('readtext.m'); %运用contour函数绘制图形readtext;for i=1:19subplot(5,4,i)contour(f(:,:,i))end运行结果如右图所示：（3）φ1=0~90，φ＝45，φ2＝0所对应的f(φ1,φ,φ2)即为f(:,10,1)。

将以下代码保存为code1_3.m文件：fopen('readtext.m');readtext;plot([0:5:90],f(:,10,1),'-bo') %运用plot函数绘制图形text(60,6,'\phi=45 \phi2=0')运行结果如下图所示：1. 编程实现以下科学计算算法，并举一例应用之。

（参考书籍《精通MATLAB 科学算法》，王正林等著，电子工业出版社，2009年）“多项式拟合”。

思考：多项式拟合是用多项式拟合曲线的一种方式，低次数下运用此方法符合较好，但较高次数下波动太大，失去真实性。

1.1 多项式曲线拟合概述对给定数据点（x i ,y i ）(i=1,2,...N ），构造m 次多项式， P （x ）=0a +mm x a x a ++ 1 （m<N ） 由曲线拟合定义，应该使得下式取极小值：∑∑=⎥⎦⎤⎢⎣⎡-Ni i ji j y x a 12通过简单的计算可得出系数是下面的线性方程组的解：⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡++m m mm m c c c c c c c c c 2112110 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡m a a a 10=⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡m b b b 10其中,c k =∑==Ni k im k x1)2,,1,0(,b k =),1,0(,m k xy k ii=∑在MATLAB 中编程实现的多项式曲线拟合函数为:multifit 功能：离散实验数据点的多项式曲线拟合。

调用格式：A=multifit(X,Y,m)其中：X 为实验数据点的x 坐标向量； Y 为实验数据点的y 坐标向量； m 为拟合多项式的次数； A 为拟合多项式的系数向量。

1.2 多项式曲线拟合编程流程图1.3 多项式曲线拟合的MATLAB 程序代码 function A=multifit(X,Y,m)%离散实验数据点的多项式曲线拟合 %实验数据点的x 坐标向量：X %实验数据点的y 坐标向量：Y%拟合多项式的次数：m输入向量X ，Y ，多项式次数m.M=N得出X 的项数为N ，Y 的项数为M建立长度为2m+1的零向量C 和长度为2m 的零向量b ，变量j ＝0,k ＝0。

j<2m+1YN显示输入不正确k<NC(j)=c(j)+X(k)^(j-1)YYj<m+2b(j)=b(j)+Y(k)*X(k)^(j-1)k=k+1建立矩阵C 将c 向量的元素依行代入将矩阵b 右除C 得到结果向量输出结果YNNN%拟合多项式的系数向量：AN=length(X);M=length(Y);if(N ~= M)disp('数据点坐标不匹配！');return;endc(1:(2*m+1))=0;b(1:(m+1))=0;for j=1:(2*m+1) %求出c和bfor k=1:Nc(j)=c(j)+X(k)^(j-1);if(j<(m+2))b(j)=b(j)+Y(k)*X(k)^(j-1); endendendC(1,:)=c(1:(m+1));for s=2:(m+1)C(s,:)=c(s:(m+s));endA=b'\C; %直接求解法求出拟合系数1.4多项式曲线拟合应用实例用二次多项式拟合下表所列的数据点。

1.4.11.4.2操作命令>> x=1:4;>> y=[4 10 18 26];>> A=multifit(x,y,2)1.4.3 输出结果输出结果为：A = 0.0489 0.1612 0.5672 即拟合的多项式为：P=0.0489+0.1612x+0.5672x21．4．4 结果如图输入向量X ，Y，多项式次数m.调用函数A=multifit(X,Y,m)输出结果２、编程解决以下科学计算问题。

2.1 问题分析解：建模：由等效电流源电路图可知各支路导纳为：Y1=1/R1+1/(j*XL); Y2=1/R2-1/(j*Xc1); Y3=1/R3-1/(j*Xc2) 均为两并联元件导纳之和，按照图中电流方向，其电流为 I1=Ua*Y1, I2=(Ub-Ua)*Y2, I3=-Ub*Y3 则a ，b 两点的电流方程为Y1Ua-Y2(Ub-Ua)=Us1/jXL+Us2/R1Y2(Ub-Ua)-Y3Ub=Us3/R3-Us4/jXc2-Us2/R2 写成矩阵形式：⎥⎦⎤⎢⎣⎡--+=⎥⎦⎤⎢⎣⎡⎥⎦⎤⎢⎣⎡+-+2/22/43/31/2/1322221R Us jXc Us R Us R Us jXL Us Ub Ua Y Y Y Y Y Y即可写成AU=B2.2 操作流程图2.3 程序代码：function fun1R1=2;R2=3;R3=4;XL=2;XC1=3;XC2=5; %给出原始数据us1=8;us2=6;us3=8;us4=15; %给出原始数据Y1=1/R1+1/( j*XL); %用复数表示各支路导路Y2=1/R2-1/( j*XC1);Y3 = 1/R3-1 /( j*XC2);A=[ Y1+ Y2,-Y2;- Y2,Y2+Y3]; %按线性方程组列ua，ub的系数矩阵B=[us1/( j*XL)+us2/R1;us3/R3+us4/(- j*XC2)-us2/R2]; %列出线性方程组右端U=A\B;ua=U(1),ub=U(2) %求ua,ubI1=ua*Y1,I2=( ub -ua)*Y2,I3=ub*Y3, %求各支路的II1R=ua/R1 ,I1L=ua/( j*XL),I2R=(ub-ua)/R2,I2C=(ub-ua)/(-j*XC1),I3R=ub/R3,I3C=ub/(- j*XC2),W=compass([ua,ub,I1,I2,I3]) %画向量图，设定此图的图柄为wset(W,'linewidth',2) %改变向量图线宽end2.4 运行结果如图：运行>> fun1ua =4.8845 - 0.5981i ub =5.4874 + 2.5752i I1 =2.1432 - 2.7413i I2 =-0.8568 + 1.2587i I3 =0.8568 + 1.7413i I1R =2.4422 - 0.2990i I1L =-0.2990 - 2.4422i I2R =0.2010 + 1.0578i I2C =-1.0578 + 0.2010iI3R =1.3718 + 0.6438iI3C =-0.5150 + 1.0975iW =179.0024180.0024181.0024182.0024183.00242.5 运行结果截图2. （2）解：由题要求，可用最小二乘拟合法拟合函数流程图程序x=[0.1 0.4 0.5 0.7 0.7 0.9];y=[0.61 0.92 0.99 1.52 1.47 2.03];cc=polyfit(x,y,2) %求出A 与B 的系数 xx=x(1):0.1:x(length(x)); yy=polyval(cc,xx); plot(xx,yy,'--') hold onplot(x,y,'x') %画出图形 axis([0,1,0,3]) xlabel('x')ylabel('y')%坐标轴名称运行结果截图。