matlab程序设计实践⽜顿法解⾮线性⽅程
中南⼤学
MATLAB程序设计实践学长有爱奉献，下载填上信息即可上交，没有下载券的⾃⾏百度。所需m⽂件照本⽂档做即可，即新建（FILE）→脚本（NEW-Sscript）→复制本⽂档代码→运⾏（会跳出保存界⾯，⽂件名默认不要修改，保存）→结果。第⼀题需要把数据⽂本⽂档和m⽂件放在⼀起。全部测试⽆误，放⼼使⽤。本⽂档针对做⽜顿法求⾮线性函数题⽬的同学，当然第⼀题都⼀样，所有⼈都可以⽤。←记得删掉这段话班级：学号：姓名：
⼀、《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 数据格式说明
解：
（1）将⽂件Data.txt 内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下：fid=fopen('data.txt'); for i=1:18tline=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;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);
将以上代码保存为readdata.m在MATLAB中运⾏，运⾏结果如下图所⽰：
将以下代码保存为code1.m⽂件：fopen('readdata.m');readdata;
[x,y,z]=meshgrid(0:5:90,0:5:90,0:5:90);slice(x,y,z,f,[45,90],[45,90],[0,45])运⾏结果如下图所⽰：
（2）将以下代码保存为code2.m⽂件：fopen('readdata.m');readdata;for i=1:19subplot(5,4,i)pcolor(f(:,:,i))nd
运⾏结果如下图所⽰：
将以下代码保存为code3.m⽂件：fopen('readdata.m');readdata;for i=1:19subplot(5,4,i)contour(f(:,:,i))end
运⾏结果如下图所⽰：
（3）φ1=0~90，φ＝45，φ2＝0所对应的f(φ1,φ,φ2)即为f(:,10,1)。将以下代码保存为code4.m⽂件：fopen('readdata.m');readdata;
plot([0:5:90],f(:,10,1),'-bo')
text(60,6,'\\phi=45 \\phi2=0')运⾏结果如下图所⽰：
⼆ 《MATLAB 程序设计实践》科学计算（24）班级： 学号： 姓名：
１、编程实现以下科学计算算法，并举⼀例应⽤之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电⼦⼯业出版社，２００９年）“⽜顿法⾮线性⽅程求解”
解：弦截法本质是⼀种割线法，它从两端向中间逐渐逼近⽅程的根；⽜顿法本质上是⼀种切线法，它从⼀端向⼀个⽅向逼近⽅程的根，其递推公式为：-=+n n x x 1)()('
n n x f x f 初始值可以取)('a f 和)('
b f 的较⼤者，这样可以加快收敛速度。 和⽜顿法有关的还有简化⽜顿法和⽜顿下⼭法。
在MATLAB 中编程实现的⽜顿法的函数为：NewtonRoot 。 功能：⽤⽜顿法求函数在某个区间上的⼀个零点。 调⽤格式：root=NewtonRoot )(```eps b a f 其中，f 为函数名；a 为区间左端点；
b 为区间右端点 eps 为根的精度；root 为求出的函数零点。 ，
⽜顿法的matlab程序代码如下：function root=NewtonRoot(f,a,b,eps)
%⽜顿法求函数f在区间[a,b]上的⼀个零点%函数名:f%区间左端点:a%区间右端点:b%根的精度:eps%求出的函数零点:rootif(nargin==3)eps=1.0e-4;end
f1=subs(sym(f),findsym(sym(f)),a);f2=subs(sym(f),findsym(sym(f)),b);if (f1==0)root=a;
endif (f2==0)root=b;endif (f1*f2>0)
disp('两端点函数值乘积⼤于0 ！');return;elsetol=1;
fun=diff(sym(f)); %求导数fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfa>dfb) %初始值取两端点导数较⼤者root=a-fa/dfa;elseroot=b-fb/dfb;end
while(tol>eps)r1=root;
fx=subs(sym(f),findsym(sym(f)),r1);
dfx=subs(sym(fun),findsym(sym(fun)),r1); %求该点的导数值 root=r1-fx/dfx; %迭代的核⼼公式tol=abs(root-r1);endend
例：求⽅程3x^2－exp（x）＝0的⼀根解：在MATLAB命令窗⼝输⼊：>> r=NewtonRoot('3*x^2-exp(x)',3,4)输出结果：X=3.7331流程图：
２、编程解决以下科学计算问题。 1)
解：这个⽅程可⽤下列步骤来解
（1）⽤eig 函数求出矩阵K-λ M 的特征值L 和特征向量U ，U 和L 满⾜
===21''00****λλU K U L IU M U
（2）在原始⽅程Mx+Kx=0两端各乘以'U 及在中间乘以对⾓矩阵'U *U ，得 'U *M*'U *U*'x +'U *K*'U *x=0
作变量置换x U x x U z z z **'21'21=??
=??=，得0*''=+z L z
这是⼀个对⾓矩阵⽅程，即可把它分两个⽅程：022''211'
'1=+=+z z z z λλ
这意味着两种振动模态可以解耦，令i i λω=2，i ω是第i 个模式的固有频率（i=1,2）。 (3)由上述的解耦模态中，给出初始条件0x ，0d x ，化为0z ，0d z ，即可求出其分量1z ，2z 。
设位置和速度的初始条件分别为[]Tx x x 02010=，[]Td d d x x x 02010=，则这三个步
骤得到的最后公式为[][][])sin 1cos .()(0021
t Mx u t Mxu u t x i d T i ii Ti
i iωωω+=∑=
系统解耦的振动模态的MATLAB 代码如下： function erziyoudu()%输⼊各原始参数
m1=input('m1=');m2=input('m2='); %质量k1=input('k1=');k2=input('k2='); %刚度%输⼊阻尼系数
流程图：
c1=input('c1=');c2=input('c2=');%给出初始条件及时间向量x0=[1;0];xd0=[0;-1];tf=50; %步数dt=0.1; %步长%构成⼆阶参数矩阵M=[m1,0;0,m2];K=[k1+k2,-k2;-k2,k2];C=[c1+c2,-c2;-c2,c2];%构成四阶参数矩阵
A=[zeros(2,2),eye(2);-M\\K,-M\\C];%四元变量的初始条件y0=[x0;xd0];
%设定计算点，作循环计算for i=1:round(tf/dt)+1t(i)=dt*(i-1);
y(:,i)=expm(A*t(i))*y0;%循环计算矩阵指数end%按两个分图绘制x1、x2曲线subplot(2,1,1),plot(t,y(1,:)),gridxlabel('t'),ylabel('y');
subplot(2,1,2),plot(t,y(2,:)),gridxlabel('t'),ylabel('y');
运⾏M⽂件，依下图所⽰在MATLAB命令窗⼝中输⼊数据：
即可得出该振动的两种模态
2)
解：第⼀步，在MATLAB命令窗⼝输⼊命令pdetool打开⼯具箱，调整x坐标范围为[0 4],y 坐标范围为[0 3].通过options选项的Axes Linits设定如下图所⽰。
第⼆步，设定矩形区域。点击⼯具箱栏中的按钮“”，拖动⿏标画出⼀矩形，并双击该矩形，设定矩形⼤⼩，如下图所⽰。
第三步，设边界条件。点击⼯具栏中的按钮“”，并双击矩形区域的相应的边线在弹出的对话框中设定边界条件。如下图所⽰，分别为各边框的边界条件。
