复合材料力学上机编程作业(计算层合板刚度)
复合材料力学上机编程作业
学院：School of Civil Engineering 专业：Engineering Mechanics
小组成员信息：James Wilson （2012031890015）、Tau Young （2012031890011）
复合材料力学学了五个星期，这是这门课的第一次编程作业。我和杨涛结成一个小组，我用的是Fortran 编制的程序，Tau Young 用的是matlab 编制。其中的算例以我的Fortran 计算结果为准。Matlab 作为可视化界面有其独到之处，在附录2中将会有所展示。
作业的内容是层合板的刚度的计算和验算，包括拉伸刚度A 、弯曲刚度D 以及耦合刚度B 。
首先要给定层合板的各个参数，具体有：层合板的层数N ；各单层的弹性常数E1、E2、21υ、G12；各单层对应的厚度；各单层对应的主方向夹角θ。然后就要计算每个单层板的二维刚度矩阵Q ，具体公式如下：
1
2
2112E E υυ=
；21
121111υυ-=
E Q ；21
12222
1υυ-=
E Q
；211211212
1υυυ-=
E Q ；
12
66G Q =
得到Q 矩阵后，根据课本上讲到的)()(11--=T Q T Q T 得到Q 。
然后根据z 坐标的定义求出0z 到n z ，接下来，最重要的一步，根据下式计算A 、B 、D 。
????
?????
-=-=-=∑∑∑=-=-=-n k k k k ij ij n k k k k ij ij n
k k k k ij ij z z Q D z z Q B z z Q A 13131
2121
1)
()(31)()(21)
()(
一、书上P110的几个问题可以归纳为以下几个类型。
（4）6层反对称角铺设层合板（T5-10）
（5）我还想验证一个书上的例题，在课本P114。三层层合板，外层厚度t1，内层10t1，正交铺设比m=0.2,。
二、验证Verchery的论文里给出的数值算例。
这里一直到Table5的数据都是从Verchery的论文中截获。
Verchery论文中的18层序列，第（21）式【laminates without bending-extension coupling】的排列有两种材料，一种是Boron-Epoxy，一种是Glass-Epoxy。而且都给出了最终的计算结果Q，A*，D*。
下面是我的Fortran计算数据文档和结果文档。
（1）Boron-Epoxy材料。
由于课本上只是分析Nx的荷载，给出了A*的数
值和课本上计算的结果一致。
结果文档
The NORMALISED extension stiffness A* equals:
2.4509E+10 4.5954E+09 -1.3893E+02
4.5954E+09 4.9017E+10 -1.2002E+03
-1.3893E+02 -1.2002E+03 8.8000E+09
The NORMALISED coupling tensor C equals:
-9.3609E+09 1.4411E-05 -4.2450E+01
1.4411E-05 9.3609E+09 -3.6673E+02
-4.2450E+01 -3.6673E+02 1.6107E-05
The NORMALISED bending stiffness D* equals:
3.3869E+10
4.5954E+09 -9.6477E+01
4.5954E+09 3.9656E+10 -8.3347E+02
-9.6477E+01 -8.3347E+02 8.8000E+09
数据文档
层数 3
层序数厚度m E1（Pa）E2（Pa）v12 v21 G12（Pa）角度（°）
1 1.00E+00 5.40E+10 1.80E+10 0.083 0.250 8.80E+09 0.00
2 1.00E+01 5.40E+10 1.80E+10 0.08
3 0.250 8.80E+09 90.00
3 1.00E+00 5.40E+10 1.80E+10 0.083 0.250 8.80E+09 0.00
（2）Glass-Epoxy材料。
（3）当然我也验证了第（22）【laminates with equal elastic properties in bending and extension】、（23）【quasi-homogeneous laminates】的排序，材料是Boron-Epoxy，下面给出计算的结果。
从下面的两个结果表中可以知道，（22）排序的确是C=0，（23）的排序的确是B=0且C=0。验证成功。
附件1：计算所用的程序代码。
PROGRAM COMPOSITE
!Coded by James Wilson
IMPLICIT NONE
REAL(8)::A(3,3),B(3,3),D(3,3),MC(5),TEMP,ROT(3,3)
!A拉伸刚度；B耦合刚度；D弯曲刚度；
!MC读入材料常数；ROT旋转矩阵
REAL(8)::TOTAL_TH,HALF_TH !总厚度；半厚度
REAL(8),ALLOCATABLE::Q(:,:,:),AL(:),T(:),Z(:),Z1(:),Z2(:),Z3(:)
!Q每层板相应刚度；AL转角；T每层板的厚度；Z坐标量INTEGER(4)::N,I,J,K,SEQ,L
!____IJK循环变量；N板的层数；SEQ序数CHARACTER(50)::CHR(8),TEMPC,filename1,filename2
!CHR、TEMPC:character variables
WRITE(*,*)"Please insert the INP file name(a.txt for example):" READ(*,*)filename1
OPEN(8,file=filename1)!Open data file
!Read in data
READ(8,*)TEMPC,N
ALLOCATE(Q(3,3,N),AL(N),T(N),Z(N+1),Z1(N),Z2(N),Z3(N)) READ(8,*)CHR(1:8)
DO I=1,N
READ(8,*)SEQ,T(I),MC(1:5),AL(I)
Q(:,:,I)=0!Calculate stiffness of each layer for the principal axis TEMP=1./(1-MC(3)*MC(4))
Q(1,1,I)=MC(1)*TEMP
Q(2,2,I)=MC(2)*TEMP
Q(3,3,I)=MC(5)
Q(1,2,I)=MC(4)*MC(2)*TEMP
Q(2,1,I)=Q(1,2,I)
AL(I)=AL(I)*3.1415926535898/180
ROT(1,1)=(cos(AL(I)))**2!Work out Rot Matrix
ROT(2,2)=ROT(1,1)
ROT(3,3)=cos(2*AL(I))
ROT(2,1)=1-ROT(1,1)
ROT(1,2)=ROT(2,1)
ROT(3,1)=0.5*sin(2*AL(I))
ROT(3,2)=-ROT(3,1)
ROT(1,3)=-2*ROT(3,1)
ROT(2,3)=-2*ROT(3,2)
Q(:,:,I)=MATMUL(MATMUL(ROT,Q(:,:,I)),TRANSPOSE(ROT)) ENDDO
TOTAL_TH=sum(T)
HALF_TH=TOTAL_TH/2
Z(1)=-HALF_TH
!Work out Z
DO I=1,N
Z(I+1)=Z(I)+T(I)
ENDDO
!calculate tensor A、B and D
DO K=1,N
Z1(K)=(Z(K+1)-Z(K))
Z2(K)=(Z(K+1)-Z(K))*(Z(K+1)+Z(K))/2
Z3(K)=(Z(K+1)**3-Z(K)**3)/3
ENDDO
A=0;B=0;D=0
WRITE(*,*)"Please insert the OUP file name(b.txt for example):" READ(*,*)filename2
OPEN(9,file=filename2)
!Write in stiffness tensor for each single ply
DO K=1,N
WRITE(9,100)K
100 FORMAT("The stiffness of number",1X,I2,2X,"ply is:")
DO I=1,3
WRITE(9,200)Q(I,:,K)
200 FORMAT(ES12.4,6X,ES12.4,6X,ES12.4)
ENDDO
WRITE(9,"(/)")
A=A+Q(:,:,K)*Z1(K)
B=B+Q(:,:,K)*Z2(K)
D=D+Q(:,:,K)*Z3(K)
ENDDO
!Output the actual stiffness of the laminate
WRITE(9,"(/)");WRITE(9,"(/)")
WRITE(9,*)"The ACTUAL stiffness tensor of the laminate:" WRITE(9,"(/)")
WRITE(9,*)"The extension stiffness A equals:"
DO I=1,3
WRITE(9,200)A(I,1:3)
ENDDO
WRITE(9,"(/)")
WRITE(9,*)"The coupling stiffness B equals:"
DO I=1,3
WRITE(9,200)B(I,1:3)
ENDDO
WRITE(9,"(/)")
WRITE(9,*)"The bending stiffness D equals:"
DO I=1,3
WRITE(9,200)D(I,1:3)
ENDDO
!Normalised tensor output
WRITE(9,"(/)");WRITE(9,"(/)")
WRITE(9,*)"The NORMALISED stiffness tensor of the laminate:" WRITE(9,"(/)")
WRITE(9,*)"The NORMALISED extension stiffness A* equals:" DO I=1,3
WRITE(9,200)A(I,1:3)/TOTAL_TH
ENDDO
WRITE(9,"(/)")
WRITE(9,*)"The NORMALISED coupling tensor C equals:"
DO I=1,3
WRITE(9,200)A(I,1:3)/TOTAL_TH-12*D(I,1:3)/TOTAL_TH**3 ENDDO
WRITE(9,"(/)")
WRITE(9,*)"The NORMALISED bending stiffness D* equals:"
DO I=1,3
WRITE(9,200)12*D(I,1:3)/TOTAL_TH**3
ENDDO
WRITE(*,*)"OUTPUT successfully,please press any key to end program!"
READ(*,*)
END PROGRAM COMPOSITE
附2
杨涛同学的MATLAB（GUI）计算程序。
主要程序：（编了个小界面，程序略长，删掉一些程序自带解释语句，添加了一些对关键语句的解释。）界面是：
作的一个算例如下：
该算例结果与组内同伴James Wilson同学基本一致，其余算例结果也基本一致，仅仅在趋近于零时有略微差异，在此不赘于。后边附上源代码：
function varargout = composit_plate(varargin)
gui_Singleton = 1;
gui_State = struct('gui_Name', mfilename, ...
'gui_Singleton', gui_Singleton, ...
'gui_OpeningFcn', @composit_plate_OpeningFcn, ...
'gui_OutputFcn', @composit_plate_OutputFcn, ...
'gui_LayoutFcn', [] , ...
'gui_Callback', []);
if nargin && ischar(varargin{1})
gui_State.gui_Callback = str2func(varargin{1});
end
if nargout
[varargout{1:nargout}] = gui_mainfcn(gui_State, varargin{:});
else
gui_mainfcn(gui_State, varargin{:});
end
function composit_plate_OpeningFcn(hObject, eventdata, handles, varargin) handles.output = hObject;
guidata(hObject, handles);
ha=axes('units','normalized','position',[0 0 1 1]);%嵌入坐标，为嵌入背景图片准备uistack(ha,'down')%作为背景
II=imread('武汉大学.jpg');%读入图片信息
image(II)
colormap hsv
set(ha,'handlevisibility','off','visible','off')
function varargout = composit_plate_OutputFcn(hObject, eventdata, handles) varargout{1} = handles.output;
function edit1_Callback(hObject, eventdata, handles)
function edit1_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit2_Callback(hObject, eventdata, handles)
function edit2_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit3_Callback(hObject, eventdata, handles)
function edit3_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit4_Callback(hObject, eventdata, ~)
function edit4_CreateFcn(hObject, eventdata, handles)
if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function pushbutton1_Callback(hObject, eventdata, handles)
syms e1e2v21g12ang%本程序采用符号运算
v12=v21*e2/e1;
q=[e1/(1-v12*v21),v21*e2/(1- v12*v21),0
v21*e2/(1-v12*v21),e2/(1-v12*v21),0
0,0,g12];
tran=[ cos(ang)^2, sin(ang)^2, -sin(2*ang)
sin(ang)^2, cos(ang)^2, sin(2*ang)
sin(2*ang)/2, -sin(2*ang)/2, cos(2*ang)];
q1=tran*q*tran';%得到Q
n=str2num(get(handles.edit2,'string'));%读入层数
n=floor(n);nn=0;
A=0;B=0;D=0;
t=str2num(get(handles.edit3,'string'));%读入每层厚度
t1=zeros(1,n+1);t1(1)=0;
for nn=1:n
t1(nn+1)=t1(nn)+t(nn);
end
ang1=str2num(get(handles.edit4,'string'));%读入每层角度
const=str2num(get(handles.edit1,'string'));%读入材料系数
t0=sum(t)/2;
t1=t1-t0;
e1=const(1);e2=const(2);v21=const(3);g12=const(4);
q11=subs(q1);
for nn=1:n
ang=ang1(nn);
nn=nn+1;
A=A+subs(q11*(t1(nn)-t1(nn-1)));
B=B+subs(0.5*q11*(t1(nn)^2-t1(nn-1)^2));
D=D+subs(1/3*q11*(t1(nn)^3-t1(nn-1)^3));
end%累加计算
set(handles.edit5,'string',num2str(A(1,:)));%以下为输出结果
set(handles.edit6,'string',num2str(A(2,:)));
set(handles.edit7,'string',num2str(A(3,:)));
set(handles.edit8,'string',num2str(B(1,:)));
set(handles.edit9,'string',num2str(B(2,:)));
set(handles.edit10,'string',num2str(B(3,:)));
set(handles.edit11,'string',num2str(D(1,:)));
set(handles.edit12,'string',num2str(D(2,:)));
set(handles.edit13,'string',num2str(D(3,:)));
function edit5_Callback(hObject, eventdata, handles) function edit5_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit6_Callback(hObject, eventdata, handles) function edit6_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit7_Callback(hObject, eventdata, handles) function edit7_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit8_Callback(hObject, eventdata, handles) function edit8_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit9_Callback(hObject, eventdata, handles) function edit9_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit10_Callback(hObject, eventdata, handles) function edit10_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit11_Callback(hObject, eventdata, handles) function edit11_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit12_Callback(hObject, eventdata, handles) function edit12_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
function edit13_Callback(hObject, eventdata, handles) function edit13_CreateFcn(hObject, eventdata, handles) if ispc && isequal(get(hObject,'BackgroundColor'),
get(0,'defaultUicontrolBackgroundColor'))
set(hObject,'BackgroundColor','white');
end
数值分析上机作业
数值分析上机实验报告 选题：曲线拟合的最小二乘法 指导老师： 专业： 学号： 姓名：
课题八曲线拟合的最小二乘法 一、问题提出 从随机的数据中找出其规律性，给出其近似表达式的问题，在生产实践和科学实验中大量存在，通常利用数据的最小二乘法求得拟合曲线。 在某冶炼过程中，根据统计数据的含碳量与时间关系，试求含碳量y 与时间t 的拟合曲线。 二、要求 1、用最小二乘法进行曲线拟合； 2、近似解析表达式为()33221t a t a t a t ++=?； 3、打印出拟合函数()t ?，并打印出()j t ?与()j t y 的误差，12,,2,1 =j ； 4、另外选取一个近似表达式，尝试拟合效果的比较； 5、*绘制出曲线拟合图*。 三、目的和意义 1、掌握曲线拟合的最小二乘法； 2、最小二乘法亦可用于解超定线代数方程组； 3、探索拟合函数的选择与拟合精度间的关系。 四、计算公式 对于给定的测量数据(x i ,f i )(i=1,2,…，n )，设函数分布为 ∑==m j j j x a x y 0)()(? 特别的，取)(x j ?为多项式 j j x x =)(? (j=0, 1,…，m )
则根据最小二乘法原理，可以构造泛函 ∑∑==-=n i m j i j j i m x a f a a a H 1 10))((),,,(? 令 0=??k a H (k=0, 1,…，m ) 则可以得到法方程 ???? ??????? ?=????????????????????????),(),(),(),(),(),(),(),(),(),(),(),(1010101111000100m m m m m m m m f f f a a a ????????????????????? 求该解方程组，则可以得到解m a a a ,,,10 ，因此可得到数据的最小二乘解 ∑=≈m j j j x a x f 0)()(? 曲线拟合：实际工作中，变量间未必都有线性关系，如服药后血药浓度与时间的关系；疾病疗效与疗程长短的关系；毒物剂量与致死率的关系等常呈曲线关系。曲线拟合是指选择适当的曲线类型来拟合观测数据，并用拟合的曲线方程分析两变量间的关系。 五、结构程序设计 在程序结构方面主要是按照顺序结构进行设计，在进行曲线的拟合时，为了进行比较，在程序设计中，直接调用了最小二乘法的拟合函数polyfit ，并且依次调用了plot 、figure 、hold on 函数进行图象的绘制，最后调用了一个绝对值函数abs 用于计算拟合函数与原有数据的误差，进行拟合效果的比较。
复合材料力学
复合材料力学 论文题目：用氧化铝填充导热和电绝缘环氧 复合材料的无缺陷石墨烯纳米片 院系班级：工程力学1302 姓名：黄义良 学号： 201314060215
用氧化铝填充导热和电绝缘环氧复合材料的无缺陷石墨烯纳米片 孙仁辉1 ，姚华1 ，张浩斌1 ，李越1 ，米耀荣2 ，于中振3 （1.北京化工大学材料科学与工程学院，有机无机复合材料国家重点实验室北京 100029;2.高级材料技术中心（CAMT ），航空航天，机械和机电工程学院J07，悉尼大学;3.北京化工大学软件物理科学与工程北京先进创新中心，北京100029） 摘要：虽然石墨烯由于其高纵横比和优异的导热性可以显着地改善聚合物的导热性，但是其导致电绝缘的严重降低，并且因此限制了其聚合物复合材料在电子和系统的热管理中的广泛应用。为了解决这个问题，电绝缘Al 2O 3用于装饰高质量（无缺陷）石墨烯纳米片（GNP ）。借助超临界二氧化碳（scCO 2），通过Al(NO 3)3 前体的快速成核和水解，然后在600℃下煅烧，在惰性GNP 表面上形成许多Al 2O 3纳米颗粒。或者，通过用缓冲溶液控制Al 2(SO 4)3 前体的成核和水解，Al 2(SO 4)3 缓慢成核并在GNP 上水解以形成氢氧化铝，然后将其转化为Al 2O 3纳米层，而不通过煅烧进行相分离。与在scCO2的帮助下的Al 2O 3@GNP 混合物相比，在缓冲溶液的帮助下制备的混合物高度有效地赋予具有优良导热性的环氧树脂，同时保持其电绝缘。具有12％质量百分比的Al 2O 3@GNP 混合物的环氧复合材料表现出1.49W /（m ·K ）的高热导率，其比纯环氧树脂高677％，表明其作为导热和电绝缘填料用于基于聚合物的功能复合材料。 关键词：聚合物复合基材料（PMCs ） 功能复合材料 电气特性 热性能 Decoration of defect-free graphene nanoplatelets with alumina for thermally conductive and electrically insulating epoxy composites Renhui Sun 1，Hua Yao 1， Hao-Bin Zhang 1，Yue Li 1，Yiu-Wing Mai 2，Zhong-Zhen Yu 3 (1.State Key Laboratory of Organic-Inorganic Composites, College of Materials Science and Engineering, Beijing University of Chemical Technology, Beijing 100029, China; 2.Centre for Advanced Materials Technology (CAMT), School of Aerospace, Mechanical and Mechatronic Engineering J07, The University of Sydney, Sydney, NSW 2006, Australia; 3.Beijing Advanced Innovation Center for Soft Matter Science and Engineering, Beijing University of Chemical Technology, Beijing 100029, China) Abstract:Although graphene can significantly improve the thermal conductivity of polymers due to its high aspect ratio and excellent thermal conductance, it causes serious reduction in electrical insulation and thus limits the wide applications of its polymer composites in the thermal management of electronics and systems. To solve this problem, electrically insulating Al 2O 3is used to decorate high quality (defect-free) graphene nanoplatelets (GNPs). Aided by supercritical carbon dioxide (scCO 2), numerous Al 2O 3 nanoparticles are formed
计算方法上机作业
计算方法上机报告 姓名： 学号： 班级： 上课班级：
说明： 本次上机实验使用的编程语言是Matlab 语言，编译环境为MATLAB 7.11.0，运行平台为Windows 7。 1. 对以下和式计算： ∑ ∞ ? ?? ??+-+-+-+=0681581482184161n n n n S n ，要求： ① 若只需保留11个有效数字，该如何进行计算； ② 若要保留30个有效数字，则又将如何进行计算； （1） 算法思想 1、根据精度要求估计所加的项数，可以使用后验误差估计，通项为： 1421114 16818485861681 n n n a n n n n n ε??= ---<< ?+++++??； 2、为了保证计算结果的准确性，写程序时，从后向前计算； 3、使用Matlab 时，可以使用以下函数控制位数： digits(位数)或vpa(变量，精度为数) （2）算法结构 1. ;0=s ?? ? ??+-+-+-+= 681581482184161n n n n t n ； 2. for 0,1,2,,n i =??? if 10m t -≤ end; 3. for ,1,2,,0n i i i =--??? ;s s t =+
（3）Matlab源程序 clear; %清除工作空间变量 clc; %清除命令窗口命令 m=input('请输入有效数字的位数m='); %输入有效数字的位数 s=0; for n=0:50 t=(1/16^n)*(4/(8*n+1)-2/(8*n+4)-1/(8*n+5)-1/(8*n+6)); if t<=10^(-m) %判断通项与精度的关系break; end end; fprintf('需要将n值加到n=%d\n',n-1); %需要将n值加到的数值 for i=n-1:-1:0 t=(1/16^i)*(4/(8*i+1)-2/(8*i+4)-1/(8*i+5)-1/(8*i+6)); s=s+t; %求和运算 end s=vpa(s,m) %控制s的精度 （4）结果与分析 当保留11位有效数字时，需要将n值加到n=7， s =3.1415926536； 当保留30位有效数字时，需要将n值加到n=22, s =3.14159265358979323846264338328。 通过上面的实验结果可以看出，通过从后往前计算，这种算法很好的保证了计算结果要求保留的准确数字位数的要求。
计算方法上机题答案
2.用下列方法求方程e^x+10x-2=0的近似根，要求误差不超过5*10的负4次方，并比较计算量 （1）二分法 （局部，大图不太看得清，故后面两小题都用局部截图） （2）迭代法
（3）牛顿法 顺序消元法 #include for(k=p+1;k 昆明理工大学工科研究生《数值分析》上机实验 学院：材料科学与工程学院 专业：材料物理与化学 学号：2011230024 姓名: 郑录 任课教师：胡杰 P277-E1 1.已知矩阵A= 10787 7565 86109 75910 ?? ?? ?? ?? ?? ??,B= 23456 44567 03678 00289 00010 ?? ?? ?? ?? ?? ?? ?? ?? ,错误！未找到引用源。 = 11/21/31/41/51/6 1/21/31/41/51/61/7 1/31/41/51/61/71/8 1/41/51/61/71/81/9 1/51/61/71/81/91/10 1/61/71/81/91/101/11?????????????????? （1）用MA TLAB函数“eig”求矩阵全部特征值。 （2）用基本QR算法求全部特征值（可用MA TLAB函数“qr”实现矩阵的QR分解）。解：MA TLAB程序如下： 求矩阵A的特征值： clear; A=[10 7 8 7;7 5 6 5;8 6 10 9;7 5 9 10]; E=eig(A) 输出结果： 求矩阵B的特征值： clear; B=[2 3 4 5 6;4 4 5 6 7;0 3 6 7 8;0 0 2 8 9;0 0 0 1 0]; E=eig(B) 输出结果： 求矩阵错误！未找到引用源。的特征值： clear; 错误！未找到引用源。=[1 1/2 1/3 1/4 1/5 1/6; 1/2 1/3 1/4 1/5 1/6 1/7; 1/3 1/4 1/5 1/6 1/7 1/8; 1/4 1/5 1/6 1/7 1/8 1/9;1/5 1/6 1/7 1/8 1/9 1/10; 1/6 1/7 1/8 1/9 1/10 1/11]; E=eig(错误！未找到引用源。) 输出结果： （2）A= 10 7877565861097 5 9 10 第一步：A0=hess(A);[Q0,R0]=qr(A0);A1=R0*Q0 返回得到： 第二部：[Q1,R1]=qr(A1);A2=R1*Q1 复合材料力学上机编程作业 学院：School of Civil Engineering专业：Engineering Mechanics 小组成员信息：James Wilson（2012031890015）、Tau Young（2012031890011）复合材料力学学了五个星期，这是这门课的第一次编程作业。我和杨涛结成一个小组，我用的是Fortran编制的程序，Tau Young用的是matlab编制。其中的算例以我的Fortran计算结果为准。Matlab作为可视化界面有其独到之处，在附录2中将会有所展示。 作业的内容是层合板的刚度的计算和验算，包括拉伸刚度A、弯曲刚度D以及耦合刚度B。 首先要给定层合板的各个参数，具体有：层合板的层数N；各单层的弹性常数 E1、E2、υ21、G12；各单层 对应的厚度；各单层对应的主方向夹角θ。然后就要计算每个单层板的二维刚度矩阵Q，具体公式如下： υ12=υ21E2 E1；Q11=E11-υ12υ21；Q22=E21-υ12υ21；Q12=υ12E1； 1-υ12υ21Q66=G12 得到Q矩阵后，根据课本上讲到的Q=(T-1)TQ(T-1)得到Q。 然后根据z坐标的定义求出z0到zn，接下来，最重要的一步，根据下式计算A、B、D。 n??Aij=∑(Qij)k(zk-zk-1) k=1??1n22?Bij=∑(Qij)k(zk-zk-1) 2k=1??1n33?Dij=∑(Qij)k(zk-zk-1)3k=1? 一、书上P110的几个问题可以归纳为以下几个类型。 第 1 页共 1 页 （4）6层反对称角铺设层合板（T5-10）第 2 页共 2 页 《数值计算方法》上机实验报告华北电力大学 实验名称数值il?算方法》上机实验课程名称数值计算方法专业班级：电力实08学生姓名：李超然学号:200801001008 成绩: 指导教师:郝育黔老师实验日期：2010年04月华北电力大学实验报告数值计算方法上机实验报吿一. 各算法的算法原理及计算机程序框图1、牛顿法求解非线性方程 *对于非线性方程，若已知根的一个近似值，将在处展开成一阶 xxfx ()0, fx ()xkk 泰勒公式 "f 0 / 2 八八,fxfxfxxxxx 0 0 0 0 0 kkkk2! 忽略高次项，有 ，fxfxfxxx 0 ()()()，，, kkk 右端是直线方程，用这个直线方程来近似非线性方程。将非线性方程的 **根代入，即fx ()0, X ，* fxfxxx 0 0 0 0, ,, kkk fx 0 fx 0 0, 解出 fX 0 *k XX,, k' fx 0 k 水将右端取为，则是比更接近于的近似值，即xxxxk, Ik, Ik fx ()k 八XX, Ikk* fx()k 这就是牛顿迭代公式。 ,2,计算机程序框图:，见, ,3,输入变量、输出变量说明: X输入变量:迭代初值，迭代精度，迭代最大次数，\0 输出变量:当前迭代次数，当前迭代值xkl ,4,具体算例及求解结果: 2/16 华北电力大学实验报吿 开始 读入 l>k /fx()0?,0 fx 0 Oxx，，01* fx ()0 XX,,，?10 kk, ,1，kN, ?xx, 10 输出迭代输出X输出奇异标志1失败标志 ,3,输入变量、输出变量说明: 结束 例：导出计算的牛顿迭代公式，并il ?算。（课本P39例2-16） 115cc （0）, 求解结果: 10. 750000 10.723837 10. 723805 10. 723805 2、列主元素消去法求解线性方程组，1,算法原理: 高斯消去法是利用现行方程组初等变换中的一种变换，即用一个不为零的数乘 -个 方程后加只另一个方程，使方程组变成同解的上三角方程组，然后再自下而上 对上三角 3/16 华北电力大学实验报告方程组求解。 列选主元是当高斯消元到第步时，从列的以下（包括）的各元素中选出绝 aakkkkkk 对值最大的，然后通过行交换将其交换到的位置上。交换系数矩阵中的 两行（包括常ekk 数项），只相当于两个方程的位置交换了，因此，列选主元不影响求解的结 ,2,计算机程序框图:，见下页, 输入变量:系数矩阵元素，常向量元素baiji 输出变量:解向量元素bbb,,12n 1、为什么结构复合材料中增强材料的形态主要为纤维？ 2、简述树脂基复合材料的优点和缺点？ 3、为什么新一代客机中复合材料用量会大幅提高？其复合材料零部件主要用到复合材料的哪些优点？ 4、为什么卫星中采用了较多的复合材料？ 答：1、利用复合材料的各种良好的力学性能用于制造结构的材料，称为结构复合材料， 它主要有基体材料和增强材料两种组分组成。其中增强材料承受主要载荷，提供复合 材料的刚度和强度，基本控制其力学性能；基体材料固定和保护增强纤维，传递纤维 间剪力和防止纤维屈曲，并改善复合材料的某些性能。用以加强制品力学性能或其他 性能的材料，在橡胶工业中又称补强剂。分纤维状和粒状材料两种。增强材料的增强 效应取决于与被增强材料的相容性，为增进相容能力，有些增强材料在使用前需要进 行表面处理。对粒状增强材料，尚需考虑其表面积（决定于粒径、形状和孔隙度）。 据报道，平均粒径在0.2μm以下的增强材料，随粒径的减小，制品的模量、抗张强度、 屈服强度和伸长率均有所增加。平均粒径较大的增强材料，由于粒径分布的不同其结 果不一致。所以，结构力学复合材料力学性能难以控制。增强材料就象树木中的纤维， 混凝土中的钢筋一样，是复合材料的重要组成部分，并起到非常重要的作用。例如在 纤维增强复合材料中，纤维是承受载荷的组元，纤维的力学性能决定了复合材料的性 能。所以说结构复合材料中增强材料的形态主要为纤维。 2、树脂基复合材料的优点：1）比强度高、比模量大2）耐疲劳性能好3）阻尼减震性 能好4）破损安全性好5）耐化学腐蚀性好6）树脂基复合材料是一种优良的电气绝缘 材料，电性能好7）树脂基复合材料热导率低、线膨胀系数小，优良的绝热材料，热 性能良好。树脂基复合材料的缺点：1)树脂基复合材料的耐热性较低2）材料的性能 分散性大。 3、用复合材料设计的飞机结构，可以推进隐身和智能结构设计的发展，有效地减少了 机体结构重量，提高了飞机运载能力，降低了发动机油耗，减少了污染排放，提高了 经济效益；复合材料优异的抗疲劳和耐介质腐蚀性能，提高了飞机结构的使用寿命和 安全性，减少了飞机的维修成本，从而提高了飞机结构的全寿命期（是指结构从论证 立项开始，有设计研制、生产研制、销售服务、使用运行、维护修理，一直到报废处 理的整个寿命期)经济性；复合材料结构有利于整个设计与整体制造技术的应用，可以 减少结构零部件的数量，提高结构的效率与可靠性，降低制造和运营成本，并可明显 改善飞机气动弹性特性，提高飞机性能。 4、正火箭导弹与航天器均要求结构重量轻,强度高。复合材料不仅兼备这两种优点,而 且还具有一些金属材料无法比拟的优良性能。卫星结构用复合材料具有重量轻、比刚 度、比强度高等特点。其碳纤维复合材料构件还具有弹性模量、热膨胀系数可设计等 特点，对卫星结构件的应用具有材料可设计的特色。 计算方法第四次上机报告 2.用欧拉方法解初值 y’=10x(1-y) 0<=x<=1 Y(0)=0 取步长h=0.1,保留5位有效数字，并与准确解相比较 分析：该题目考察欧拉方法解初值问题 程序如下： function Heun(a,b,y0,n) h=(b-a)/n;x=a:h:b; y=y0*ones(1,n+1); for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1)); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc); end for k=1:n+1 fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end function z=f(xx,yy) z=10*xx*(1-yy); 运行结果： >> Heun(0,1,0,10) x[0]=0.000000 y[0]=0.000000 x[1]=0.100000 y[1]=0.050000 x[2]=0.200000 y[2]=0.183000 x[3]=0.300000 y[3]=0.362740 x[4]=0.400000 y[4]=0.547545 x[5]=0.500000 y[5]=0.705905 x[6]=0.600000 y[6]=0.823543 x[7]=0.700000 y[7]=0.901184 x[8]=0.800000 y[8]=0.947627 x[9]=0.900000 y[9]=0.973290 x[10]=1.000000 y[10]=0.986645 >> 分析： 该结果与准确结果相比比较接近，但是有一定的误差。 6.用四阶龙格—库塔公式解第三题中的初值问题，取步长h=0.2，保留五位有效数字。 题目目的分析： 该题考查四阶龙格-库塔方法和改进欧拉方法求解精确度问题。 程序： 改进欧拉法： function Heun(a,b,y0,n) h=(b-a)/n;x=a:h:b; y=y0*ones(1,n+1); for j=2:n+1 yp=y(j-1)+h*f(x(j-1),y(j-1)); yc=y(j-1)+h*f(x(j),yp); y(j)=1/2*(yp+yc); end for k=1:n+1 fprintf('x[%d]=%f\ty[%d]=%f\n',k-1,x(k),k-1,y(k)); end 计算方法 一、填空题 1.假定x ≤1，用泰勒多项式?+??+++=! !212n x x x e n x ,计算e x 的值，若要求截断误差不超过0.005，则n=_5___ 2. 解 方 程 03432 3=-+x － x x 的牛顿迭代公式 )463/()343(121121311+--+--=------k k k k k k k x x x x x x x 3.一阶常微分方程初值问题 ?????= ='y x y y x f y 0 0)() ,(,其改进的欧拉方法格式为)],(),([21 1 1 y x y x y y i i i i i i f f h +++++= 4.解三对角线方程组的计算方法称为追赶法或回代法 5. 数值求解初值问题的四阶龙格——库塔公式的局部截断误差为o(h 5 ) 6.在ALGOL 中，简单算术表达式y x 3 + 的写法为x+y ↑3 7.循环语句分为离散型循环,步长型循环,当型循环. 8.函数)(x f 在[a,b]上的一次（线性）插值函数= )(x l )()(b f a b a x a f b a b x --+-- 9.在实际进行插值时插值时，将插值范围分为若干段，然后在每个分段上使用低阶插值————如线性插值和抛物插值，这就是所谓分段插值法 10、数值计算中，误差主要来源于模型误差、观测误差、截断误差和舍入误差。 11、电子计算机的结构大体上可分为输入设备 、 存储器、运算器、控制器、 输出设备 五个主要部分。 12、算式2 cos sin 2x x x +在ALGOL 中写为))2cos()(sin(2↑+↑x x x 。 13、ALGOL 算法语言的基本符号分为 字母 、 数字 、 逻辑值、 定义符四大 东南大学数值分析上机作业 汇总 -标准化文件发布号：（9456-EUATWK-MWUB-WUNN-INNUL-DDQTY-KII 数值分析上机报告 院系： 学号： 姓名： 目录 作业1、舍入误差与有效数 (1) 1、函数文件cxdd.m (1) 2、函数文件cddx.m (1) 3、两种方法有效位数对比 (1) 4、心得 (2) 作业2、Newton迭代法 (2) 1、通用程序函数文件 (3) 2、局部收敛性 (4) （1）最大δ值文件 (4) （2）验证局部收敛性 (4) 3、心得 (6) 作业3、列主元素Gauss消去法 (7) 1、列主元Gauss消去法的通用程序 (7) 2、解题中线性方程组 (7) 3、心得 (9) 作业4、三次样条插值函数 (10) 1、第一型三次样条插值函数通用程序： (10) 2、数据输入及计算结果 (12) 作业1、舍入误差与有效数 设∑ =-=N j N j S 2 2 11 ，其精确值为?? ? ??---1112321N N . （1）编制按从小到大的顺序1 1 131121222-? ??+-+-=N S N ，计算N S 的通用程序； （2）编制按从大到小的顺序()1 21 11111222-???+--+-=N N S N ，计算N S 的通用程序； （3）按两种顺序分别计算642101010,,S S S ,并指出有效位数； （4）通过本上机你明白了什么？ 程序： 1、函数文件cxdd.m function S=cxdd(N) S=0; i=2.0; while (i<=N) S=S+1.0/(i*i-1); i=i+1; end script 运行结果（省略>>）： S=cxdd(80) S= 0.737577 2、函数文件cddx.m function S=cddx (N) S=0; for i=N:-1:2 S=S+1/(i*i-1); end script 运行结果（省略>>）： S=cddx(80) S= 0.737577 3、两种方法有效位数对比 复合材料力学上机作业 （2013年秋季） 班级力学C102 学生姓名赵玉鹰 学号105634 成绩 河北工业大学机械学院 2013年12月30日 作业1 单向板刚度及柔度的计算 一、要 求 （1）选用FORTRAN 、VB 、MAPLE 或MATLAB 编程计算下列各题； （2）上机报告内容：源程序、题目内容及计算结果； （3）材料工程常数的数值参考教材自己选择； （4）上机学时：2学时。 二、题 目 1、已知单层板材料工程常数1E ，2E ，12G ，计算柔度矩阵[S ]和刚度矩阵[Q ]。（玻璃/环氧树脂单层板材料的MPa 1090.341?=E ，MPa 1030.142?=E ，MPa 1042.0412?=G ，25.021=μ，MPa 1001=σ，MPa 302-=σ，MPa 1012=τ） ●Maple 程序 > restart: > with(linalg): > E[1]:=3.9e10: > E[2]:=1.3e10: > G[12]:=0.42e10: > mu[21]:=0.25: > mu[12]:=E[1]*mu[21]/E[2]: > Q[11]:=E[1]/(1-mu[12]*mu[21]): > Q[12]:=mu[12]*E[2]/(1-mu[12]*mu[21]): > Q[13]:=0: > Q[21]:=Q[12]: > Q[22]:=E[2]/(1-mu[12]*mu[21]): > Q[23]:=0: > Q[31]:=Q[13]: > Q[32]:=Q[23]: > Q[33]:=G[12]: >Q:=evalf(matrix(3,3,[[Q[11],Q[12],Q[13]],[Q[21],Q[22], Q[23]],[Q[31],Q[32],Q[33]]]),4); 计算方法实验报告 班级： 学号： 姓名： 成绩： 1 舍入误差及稳定性 一、实验目的 （1）通过上机编程，复习巩固以前所学程序设计语言及上机操作指令； （2）通过上机计算，了解舍入误差所引起的数值不稳定性 二、实验内容 1、用两种不同的顺序计算10000 21n n -=∑，分析其误差的变化 2、已知连分数() 1 01223//(.../)n n a f b b a b a a b =+ +++，利用下面的算法计算f ： 1 1 ,i n n i i i a d b d b d ++==+ (1,2,...,0 i n n =-- 0f d = 写一程序，读入011,,,...,,,...,,n n n b b b a a 计算并打印f 3、给出一个有效的算法和一个无效的算法计算积分 1 041 n n x y dx x =+? (0,1,...,1 n = 4、设2 2 11N N j S j == -∑ ，已知其精确值为1311221N N ?? -- ?+?? （1）编制按从大到小的顺序计算N S 的程序 （2）编制按从小到大的顺序计算N S 的程序 （3）按两种顺序分别计算10001000030000,,,S S S 并指出有效位数 三、实验步骤、程序设计、实验结果及分析 1、用两种不同的顺序计算10000 2 1n n -=∑，分析其误差的变化 （1）实验步骤： 分别从1~10000和从10000~1两种顺序进行计算，应包含的头文件有stdio.h 和math.h （2）程序设计： a.顺序计算 #include 1． 设?+=1 05dx x x I n n ， （1） 由递推公式n I I n n 1 51+-=-，从0I 的几个近似值出发，计算20I ； 解：易得：0I =ln6-ln5=0.1823, 程序为： I=0.182; for n=1:20 I=(-5)*I+1/n; end I 输出结果为：20I = -3.0666e+010 （2） 粗糙估计20I ，用n I I n n 51 5111+- =--，计算0I ； 因为 0095.05 6 0079.01020 201 020 ≈<<≈??dx x I dx x 所以取0087.0)0095.00079.0(2 1 20=+= I 程序为：I=0.0087; for n=1:20 I=(-1/5)*I+1/(5*n); end I 0I = 0.0083 （3） 分析结果的可靠性及产生此现象的原因(重点分析原因)。 首先分析两种递推式的误差；设第一递推式中开始时的误差为000I I E '-=，递推过程的舍入误差不计。并记n n n I I E '-=，则有01)5(5E E E n n n -==-=-Λ。因为=20E 20020)5(I E >>-，所此递推式不可靠。而在第二种递推式中n n E E E )5 1(5110-==-=Λ，误差在缩小， 所以此递推式是可靠的。出现以上运行结果的主要原因是在构造递推式过程中，考虑误差是否得到控制， 即算法是否数值稳定。 2． 求方程0210=-+x e x 的近似根，要求4 1105-+?<-k k x x ，并比较计算量。 （1） 在[0，1]上用二分法； 程序：a=0;b=1.0; while abs(b-a)>5*1e-4 c=(b+a)/2; 《复合材料力学》课程上机指导书（力学101，力学C101-2） 河北工业大学机械学院力学系 2013年9月 目录 作业1单向板刚度及柔度的计算(2学时) (1) 作业2单向板的应力、应变计算(2学时) (2) 作业3绘制表观工程常数随 的变化规律(3学时) (3) 作业4绘制强度准则的理论曲线(包络线)(3学时) (4) 作业5层合板的刚度计算(3学时) (5) *作业6层合板的强度计算(4学时) (6) 附录作业提交说明……………………………………………. . 7 注：带“*”的题目可作为自愿选做题。 作业1 单向板刚度及柔度的计算 一、要 求 （1）选用FORTRAN 、VB 、MAPLE 或MATLAB 编程计算下列各题； （2）上机报告内容：源程序、题目内容及计算结果； （3）材料工程常数的数值参考教材自己选择； （4）上机学时：2学时。 二、题 目 1、已知单层板材料工程常数1E ，2E ，12G ，计算柔度矩阵[S ]和刚度矩阵[Q ]。（玻璃/环氧树脂单层板材料的MPa 1090.341?=E ，MPa 1030.142?=E ， MPa 1042.0412?=G ，25.021=μ， MPa 1001=σ，MPa 302-=σ，MPa 1012=τ） 2、已知单层板材料工程常数1E ，2E ，12G ，21μ及θ，计算柔度矩阵][S 和刚度矩阵][Q 。（M P a 1090.341?=E ，MPa 1030.142?=E ，MPa 1042.0412?=G ，25.021=μ，?=30θ） 作业2 单向板的应力、应变计算 一、要 求 1、选用FORTRAN 、VB 、MAPLE 或MATLAB 编程计算下列各题； 2、上机报告内容：源程序、题目内容及计算结果； 3、材料工程常数的数值参考教材自己选择； 4、上机学时：2学时。 二、题 目 1、已知单向板的应力x σ、y σ、xy τ，工程常数1E ，2E ，12G ，21μ及θ，求x ε、 y ε、xy γ；1σ、2σ、12τ；1ε、2ε、12γ。 （知?=30θ，应力MPa 160=x σ，MPa 60=y σ，MPa 20=xy τ，工程常数MPa 1090.341?=E ，MPa 1030.142?=E ，MPa 1042.0412?=G ，25.021=μ，?=30θ） 2、已知1σ、2σ、12τ，工程常数1E ，2E ，12G ，21μ及θ，求1ε、2ε、12γ；x ε、y ε、 xy γ；x σ、y σ、xy τ。 （知MPa 1001=σ，MPa 302-=σ，MPa 1012=τ，MPa 1090.341?=E ，MPa 1030.142?=E ，MPa 1042.0412?=G ，25.021=μ，?=30θ） 2015.1.9 上机作业题报告 JONMMX 2000 1.Chapter 1 1.1题目 设S N =∑1j 2?1 N j=2 ，其精确值为 )1 1 123(21+--N N 。 （1）编制按从大到小的顺序1 1 131121222-+ ??+-+-=N S N ，计算S N 的通用程序。 （2）编制按从小到大的顺序1 21 1)1(111222-+ ??+--+-= N N S N ，计算S N 的通用程序。 （3）按两种顺序分别计算64210,10,10S S S ,并指出有效位数。（编制程序时用单精度） （4）通过本次上机题，你明白了什么？ 1.2程序 1.3运行结果 1.4结果分析 按从大到小的顺序，有效位数分别为：6，4，3。 按从小到大的顺序，有效位数分别为：5，6，6。 可以看出，不同的算法造成的误差限是不同的，好的算法可以让结果更加精确。当采用从大到小的顺序累加的算法时，误差限随着N 的增大而增大，可见在累加的过程中，误差在放大，造成结果的误差较大。因此，采取从小到大的顺序累加得到的结果更加精确。 2.Chapter 2 2.1题目 （1）给定初值0x 及容许误差ε，编制牛顿法解方程f(x)=0的通用程序。 （2）给定方程03 )(3 =-=x x x f ,易知其有三个根3,0,3321= *=*-=*x x x ○1由牛顿方法的局部收敛性可知存在,0>δ当),(0δδ+-∈x 时，Newton 迭代序列收敛于根x2*。试确定尽可能大的δ。 ○2试取若干初始值，观察当),1(),1,(),,(),,1(),1,(0+∞+-----∞∈δδδδx 时Newton 序列的收敛性以及收敛于哪一个根。 （3）通过本上机题，你明白了什么？ 2.2程序 第一次&第二次上机作业 上机作业： 1.在Matlab上执行：>> 5.1-5-0.1和>> 1.5-1-0.5 给出执行结果，并简要分析一下产生现象的原因。 解：执行结果如下： 在Matlab中，小数值很难用二进制进行描述。由于计算精度的影响，相近两数相减会出现误差。 2.（课本181页第一题） 解：（1）n=0时，积分得I0=ln6-ln5，编写如下图代码 从以上代码显示的结果可以看出，I 20的近似值为0.7465 （2）I I =∫I I 5+I 10dx,可得∫I I 610dx ≤∫I I 5+I 10dx ≤∫I I 510dx,得 16(I +1)≤I I ≤15(I +1),则有1126≤I 20≤1105, 取I 20=1 105 ,以此逆序估算I 0。代码段及结果如下图所示 （3）从I20估计的过程更为可靠。首先根据积分得表达式是可知，被积函数随着n的增大，其所围面积应当是逐步减小的，即积分值应是随着n的递增二单调减小的，（1）中输出的值不满足这一条件，（2）满足。设I I表示I I的近似值，I I-I I=(?5)I(I0?I0)(根据递推公式可以导出此式),可以看出，随着n的增大，误差也在增大，所以顺序估计时，算法不稳定性逐渐增大，逆序估计情况则刚好相反，误差不断减小，算法逐渐趋于稳定。 2.（课本181页第二题） （1）上机代码如图所示 求得近似根为0.09058 （2）上机代码如图所示 得近似根为0.09064； （3）牛顿法上机代码如下 计算所得近似解为0.09091 第三次上机作业上机作业181页第四题 线性方程组为 [1.13483.8326 0.53011.7875 1.16513.4017 2.53301.5435 3.4129 4.9317 1.23714.9998 8.76431.3142 10.67210.0147 ][ I1 I2 I3 I4 ]=[ 9.5342 6.3941 18.4231 16.9237 ] （1）顺序消元法 A=[1.1348,3.8326,1.1651,3.4017;0.5301,1.7875,2.5330,1.5435; 3.4129, 4.9317,8.7643,1.3142;1.2371,4.9998,10.6721,0.0147]; b=[9.5342;6.3941;18.4231;16.9237]; 上机代码（函数部分）如下 function [b] = gaus( A,b )%用b返回方程组的解 B=[A,b]; n=length(b); RA=rank(A); RB=rank(B); 数值计算方法上机题目1 1、实验1. 病态问题 实验目的： 算法有“优”与“劣”之分，问题也有“好”和“坏”之别。所谓坏问题就是问题本身的解对数据变化的比较敏感，反之属于好问题。希望读者通过本实验对此有一个初步的体会。 数值分析的大部分研究课题中，如线性代数方程组、矩阵特征值问题、非线性方程及方程组等都存在病态的问题。病态问题要通过研究和构造特殊的算法来解决，当然一般要付出一些代价（如耗用更多的机器时间、占用更多的存储空间等）。 问题提出： 考虑一个高次的代数多项式 ∏=-= ---=20 1)()20)...(2)(1()(k k x x x x x p （E1-1） 显然该多项式的全部根为l ，2，…，20，共计20个，且每个根都是单重的（也称为简 单的）。现考虑该多项式方程的一个扰动 0)(19 =+x x p ε (E1-2) 其中ε是一个非常小的数。这相当于是对（E1-1）中19 x 的系数作一个小的扰动。我们希望比较（E1-1）和（E1-2）根的差别，从而分析方程（E1-1）的解对扰动的敏感性。 实验内容： 为了实现方便，我们先介绍两个 Matlab 函数：“roots ”和“poly ”，输入函数 u ＝roots （a ） 其中若变量a 存储1+n 维的向量，则该函数的输出u 为一个n 维的向量。设a 的元素依次为121,...,,+n a a a ，则输出u 的各分量是多项式方程 0...1121=++++-n n n n a x a x a x a 的全部根，而函数 b=poly(v) 的输出b 是一个n ＋1维变量，它是以n 维变量v 的各分量为根的多项式的系数。可见“roots ”和“Poly ”是两个互逆的运算函数. ve=zeros(1,21); ve(2)=ess; roots(poly(1:20))+ve) 上述简单的Matlab 程序便得到（E1-2）的全部根，程序中的“ess ”即是（E1-2）中的ε。 实验要求： （1）选择充分小的ess ，反复进行上述实验，记录结果的变化并分析它们。如果扰动项的系数ε很小，我们自然感觉（E1-1）和(E1-2)的解应当相差很小。计算中你有什么出乎意料的发现？表明有些解关于如此的扰动敏感性如何？ （2）将方程（E1-2）中的扰动项改成18 x ε或其他形式，实验中又有怎样的现象出现？数值分析上机作业
复合材料力学上机编程作业(计算层合板刚度)要点
《数值计算方法》上机实验报告
复合材料力学设计作业1
计算方法上机作业
计算方法试题库讲解
东南大学数值分析上机作业汇总
复合材料力学大作业
计算方法上机实习题大作业(实验报告).
(完整版)数值计算方法上机实习题答案
2013秋复合材料力学上机作业
东南大学-数值分析上机题作业-MATLAB版
计算方法上机作业集合
数值分析上机作业1-1