中南大学材料专业MATLAB实践4号分析.docx
- 文档编号:27161008
- 上传时间:2023-06-27
- 格式:DOCX
- 页数:24
- 大小:1.50MB
中南大学材料专业MATLAB实践4号分析.docx
《中南大学材料专业MATLAB实践4号分析.docx》由会员分享,可在线阅读,更多相关《中南大学材料专业MATLAB实践4号分析.docx(24页珍藏版)》请在冰豆网上搜索。
中南大学材料专业MATLAB实践4号分析
一、《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分布情况。
备注:
data.txt数据格式说明
将文件Data.txt内的数据按照要求读取到矩阵f(phi1,phi,phi2)中,
代码如下:
fid=fopen('data.txt');%读取数据文件Data.txt
fori=1:
18
tline=fgetl(fid);
end
phi1=1;phi=1;phi2=1;line=0;
f=zeros(19,19,19);
while~feof(fid)
tline=fgetl(fid);
data=str2num(tline);
line=line+1;
ifmod(line,20)==1
phi2=(data/5)+1;
phi=1;
else
forphi1=1:
19
f(phi1,phi,phi2)=data(phi1);
end
phi=phi+1;
end
end
fclose(fid);
建立新的脚本,将以上代码保存为m文件,命名为readtext.m,在MATLAB运行后。
其变量f结果为:
(1).用Slice函数给出其整体分布特征
程序代码:
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函数绘制图形
将程序代码保存为code1_1.m文件,其运行结果:
(2).
.用subplot函数和pcolor函数给出(φ2=0,5,10,15,20,25,30,35…90)切面上f分布情况;
程序代码:
fopen('readtext.m');
readtext;fori=1:
19
subplot(5,4,i)
pcolor(f(:
:
i))%运用pcolor函数绘制图形
end
将上述代码保存为code1_2_1.m文件,其运行结果:
.使用subplot函数和contour函数给出(φ2=0,5,10,15,20,25,30,35…90)切面上f分布情况;
程序代码:
fopen('readtext.m');%运用contour函数绘制图形readtext;
fori=1:
19
subplot(5,4,i)
contour(f(:
:
i))
end
将上述代码保存为code1_2_2.m文件,运行结果:
(3).用plot函数给出沿α取向线(φ1=0~90,φ=45,φ2=0)的f分布情况。
程序代码:
fopen('readtext.m');
readtext;
plot([0:
5:
90],f(:
10,1),'-bo')%运用plot函数绘制图形
text(60,6,'\phi=45\phi2=0')
将上述代码保存为code1_3.m文件,运行结果:
二《MATLAB程序设计实践》科学计算(04)
班级
1、编程实现以下科学计算算法,并举一例应用之。
(参考书籍《精通MALAB科学计算》,王正林等著,电子工业出版社,2009年)“切比雪夫逼近”
算法说明:
当一个连续函数定义在区间[-1,1]上时,它可以展开成切比雪夫级数。
即:
其中
为
次切比雪夫多项式,具体表达可通过递推得出:
,
它们之间满足如下的正交关系:
在实际应用中,可根据所需的精度来截取有限的项数,切比雪夫级数中的系数由下式决定:
流程图:
程序源代码(m文件):
functionf=Chebyshev(y,k,x0)
%用切比雪夫多项式逼近已知函数
%已知函数:
y
%逼近已知函数所需项数:
k
%逼近点的x坐标:
x0
%求得的切比雪夫逼近多项式或在x0处的逼近值:
f
symst;%定义sym型变量t
T(1:
k+1)=t;
T
(1)=sym('1');
T
(2)=t;%定义切比雪夫多项式矩阵T(n),并规定前两项为1,t
c(1:
k+1)=sym('0');
c
(1)=int(subs(y,findsym(sym(y)),sym('t'))*T
(1)/sqrt(1-t^2),t,-1,1)/pi;
c
(2)=2*int(subs(y,findsym(sym(y)),sym('t'))*T
(2)/sqrt(1-t^2),t,-1,1)/pi;
%根据前两项多项式计算其系数c
(1),c
(2)
f=c
(1)+c
(2)*t;%由第一项和第二项多项式及其系数计算切比雪夫逼近值
fori=3:
k+1%从第3项到第k+1项进行循环
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-t^2),t,-1,1)/2;
%计算第i项系数
f=f+c(i)*T(i);%计算切比雪夫逼近值
f=vpa(f,6);
if(i==k+1)
if(nargin==3)
f=subs(f,'t',x0);%输入三项,则输出x0处六级精度的切比雪夫逼近值
end
f=vpa(f,6);
end
end
根据代码可知,在MATLAB中编程实现的切比雪夫逼近法函数为:
Chebyshev。
功能:
用切比雪夫多项式逼近已知函数。
调用格式:
其中,y为已知函数;
k为逼近已知函数所需项数;
f是求得的切比雪夫逼近多项式或在x0处的逼近值。
应用实例:
切比雪夫应用实例。
用切比雪夫公式(取6项)逼近函数
,并求当x=1时的函数值。
解:
利用程序求解方程,在MATLAB命令窗口中输入:
>>Chebyshev('1/(2-x)',5)%调用创建的函数,输出切比雪夫多项式的5个项
结果:
ans=
0.277013*t+0.0647768*t*(2.0*t^2-1.0)-0.00501051*t*(2.0*t*(t-2.0*t*(2.0*t^2-1.0))+2.0*t^2-1.0)-0.0186995*t*(t-2.0*t*(2.0*t^2-1.0))+0.24175*t^2+0.456475
再在MATLAB命令窗口中输入:
>>Chebyshev('1/(2-x)',5,1)%调用创建的函数,输出当x=1时的函数值逼近值
结果:
ans=
1.06372
在MATLAB命令窗口计算符号表达式x=1的极限:
>>limit(sym('1/(2-x)'),'x',1)%计算其真值,并与其进行比较
结果:
ans=
1
输出结果:
绝对误差:
1.06372-1=0.06372;相对误差:
0.06372/1*100%=6.372%
2、编程解决以下科学计算问题。
1)
问题分析:
1.相关参数字母设定:
如图,由在P的作用,点A移动到点A’,设在水平方向移动距离为dx,在竖直方向移动距离为dy,P与水平的夹角为α。
设各杆的轴力为Ni(i=1,2,3,…n);各杆与水平面夹角分别αi(i=1,2,3,…n);各杆的变化长度分别为dLi(i=1,2,3,…n);
2.解题思路:
(1)取一根杆分析,由题意可知每根杆的伸长量为dx和dy在杆上的分量;故可得方程
:
(2)将各轴力和外力P水平(x)竖直(y)正交分解,由受力平衡可得:
x方向有方程
:
y方向有方程
:
(3)联立上述
可得n+2等式,可解得n个轴力,以及A点位移dx和dy,联立的线性方程组如下:
.
.
.
(4)建立[Pcosα,Psinα,0,0,0…0]’的常数矩阵,以及如下所示的系数矩阵:
(5)再用求逆法求解此线性方程组,即用常数矩阵除以系数矩阵,得出结果。
程序源代码:
Ei=input('请输入各杆的刚度:
(注意用[]括起来)');%输入刚度矩阵Ei
Li=input('请输入各杆的长度:
(注意用[]括起来)');%输入杆的长度矩阵Li
Ai=input('请输入各杆的横截面积:
(注意用[]括起来)');%输入杆的横截面积矩阵Ai
ai=input('请输入各杆与水平面的夹角:
(注意用[]括起来)');%输入杆与水平面的夹角矩阵ai
P=input('请输入外力P:
');%输入外加力P
a=input('请输入P与水平面的夹角:
');%输入外加力P与x的夹角
n1=length(Ei);n2=length(Li);n3=length(Ai);n4=length(ai);
if(n1~=n2|n1~=n3|n1~=n4)
disp('输入数据错误')
else
n=n1;
end%判断赋值维度是否一致
Ki=Li./(Ei.*Ai);
C=zeros(n+2,1);
C(1,1)=P*cos(a);
C(2,1)=P*sin(a);
C(3:
n+2,1)=zeros(n,1);%建立方程组等号右边常数的矩阵
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);
end
D(3:
n+2,n+1)=(-cos(ai));
D(3:
n+2,n+2)=(-sin(ai));%建立方程组系数矩阵
x=D\C;x=x';%求解该线性方程组,得出位移以及每根杆的轴力
disp('节点在x、y方向上的位移分别:
')
x(n+1:
n+2)
disp('各杆的轴力分别为:
')
x(1:
n)%输出结果
流程图:
应用实例:
设三根杆组成的支架如下图所示,挂一重物P=3000N。
设L=3m,各杆的横截面积分别为:
A1=150
10-6m2,A2=200
10-6m2,A3=300
10-6m2,材料的弹性模量均为E=200
109N/m2,求各杆所受力的大小以及C点位移
解:
将前面的程序源代码保存为main.m文件运行,
根据提示依次输入题中数据:
[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;
输出结果:
节点在x、y方向上的位移分别:
ans=
1.0e-03*
0.33990.0420
各杆的轴力分别为:
ans=
1.0e+03*
1.78650.5595-2.9794
2)
流程图:
(1)
程序代码:
>>fun=inline('(1./((2*pi).^0.5)).*exp(-x.^2./2)');
>>quad(fun,0,1)
在MATLAB的命令行窗口输入上述代码,结果:
ans=
0.3413
(2)
程序代码:
>>symsx
>>a=int(sin(x)/x,x,0,1)
>>vpa(a,5)
在MATLAB的命令行窗口输入上述代码,结果:
ans=
0.94608
(3)
;
程序代码:
>>fun=inline('x.^-x');
>>quad(fun,0,1)
在MATLAB的命令行窗口输入上述代码,结果:
ans=
1.2913
(4)
;
解:
程序代码:
>>int(sym('exp(2*x)*[sin(x)]^2'),0,2*pi)
在MATLAB的命令行窗口输入上述代码,结果:
ans=
exp(4*pi)/8-1/8
(5)
;
解:
将r用x代替,θ用y代替;
则程序代码:
>>fun=inline('(1+(x.^2).*sin(y)).^0.5');
>>dblquad(fun,0,1,0,2*pi)
在MATLAB的命令行窗口输入上述代码,结果:
ans=
6.1879
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 中南 大学 材料 专业 MATLAB 实践 分析