matlab程序设计实践-牛顿法解非线性方程

tol=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);
备注：数据格式说明
解：
（1）将文件内的数据按照要求读取到矩阵f(phi1,phi,phi2)中，代码如下：
fid=fopen('');
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
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)+1
t(i)=dt*(i-1);
y(:,i)=expm(A*t(i))*y0;%循环计算矩阵指数
班级：
学号：
姓名：
一、《LAB程序设计实践》Matlab基础
表示多晶体材料织构的三维取向分布函数（f＝f（φ1，φ，φ2））是一个非常复杂的函数，难以精确的用解析函数表达，通常采用离散空间函数值来表示取向分布函数，是三维取向分布函数的一个实例。由于数据量非常大，不便于分析，需要借助图形来分析。请你编写一个matlab程序画出如下的几种图形来分析其取向分布特征：
root为求出的函数零点。
，
牛顿法的matlab程序代码如下：
functionroot=NewtonRoot(f,a,b,eps)
%牛顿法求函数f在区间[a,b]上的一个零点
%函数名:f
%区间左端点:a
%区间右端点:b
%根的精度:eps
%求出的函数零点:root
if(nargin==3)
eps=;
end
初始值可以取 和 的较大者，这样可以加快收敛速度。
和牛顿法有关的还有简化牛顿法和牛顿下山法。
在MATLAB中编程实现的牛顿法的函数为：NewtonRoot。
功能：用牛顿法求函数在某个区间上的一个零点。
调用格式：root=NewtonRoot
其中， 为函数名；
为区间左端点；
为区间右端点
eps为根的精度；
end
%按两个分图绘制x1、x2曲线
subplot(2,1,1),plot(t,y(1,:)),grid
xlabel('t'),ylabel('y');
subplot(2,1,2),plot(t,y(2,:)),grid
xlabel('t'),ylabel('y');
运行M文件，依下图所示在MATLAB命令窗口中输入数据：
f1=subs(sym(f),findsym(sym(f)),a);
f2=subs(sym(f),findsym(sym(f)),b);
if(f1==0)
root=a;
end
if(f2==0)
root=b;
end
if(f1*f2>0)
disp('两端点函数值乘积大于0 ！');
return;
else
二 《MATLAB程序设计实践》科学计算（24）
班级：
学号：
姓名：
１、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）
“牛顿法非线性方程求解”
解：弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，它从一端向一个方向逼近方程的根，其递推公式为：
（2）将以下代码保存为文件：
fopen('');
readdata;
fori=1:19
subplot(5,4,i)
pcolor(f(:,:,i))
nd
运行结果如下图所示：
将以下代码保存为文件：
fopen('');
readdata;
fori=1:19
subplot(5,4,i)
contour(f(:,:,i))
即可得出该振动的两种模态
2)
解：第一步，在MATLAB命令窗口输入命令pdetool打开工具箱，调整x坐标范围为[04],y坐标范围为[03].通过options选项的Axes Linits设定如下图所示。
第二步，设定矩形区域。点击工具箱栏中的按钮“ ”，拖动鼠标画出一矩形，并双击该矩形，设定矩形大小，如下图所示。
中南大学
MATLAB程序设计实践
学长有爱奉献LE）→脚本（NEW-Sscript）→复制本文档代码→运行（会跳出保存界面，文件名默认不要修改，保存）→结果。第一题需要把数据文本文档和m文件放在一起。全部测试无误，放心使用。本文档针对做牛顿法求非线性函数题目的同学，当然第一题都一样，所有人都可以用。←记得删掉这段话
第三步，设边界条件。点击工具栏中的按钮“ ”，并双击矩形区域的相应的边线在弹出的对话框中设定边界条件。如下图所示，分别为各边框的边界条件。
第四步，设定方程。单击工具栏中的按钮“ ”，在PDE模式下选择方程类型，如下图所示，并在其中设定参数。
第五步，单击工具栏中的按钮“ ”，拆分区域为若干子区域，如下图所示。
（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分布情况。
k1=input('k1=');k2=input('k2='); %刚度
%输入阻尼系数
c1=input('c1=');c2=input('c2=');
%给出初始条件及时间向量
x0=[1;0];
xd0=[0;-1];
tf=50; %步数
dt=; %步长
%构成二阶参数矩阵
M=[m1,0;0,m2];
K=[k1+k2,-k2;-k2,k2];
(3)由上述的解耦模态中，给出初始条件 ， ，化为 ， ，即可求出其分量 ， 。
设位置和速度的初始条件分别为 ， ，则这三个步 骤得到的最后公式为
系统解耦的振动模态的MATLAB代码如下：
function erziyoudu()
%输入各原始参数
m1=input('m1=');m2=input('m2='); %质量
end
运行结果如下图所示：
（3）φ1=0~90，φ＝45，φ2＝0所对应的f(φ1,φ,φ2)即为f(:,10,1)。将以下代码保存为文件：
fopen('');
readdata;
plot([0:5:90],f(:,10,1),'-bo')
text(60,6,'\phi=45 \phi2=0')
运行结果如下图所示：
root=r1-fx/dfx;%迭代的核心公式
tol=abs(root-r1);
end
end
例：求方程3x^2－exp（x）＝0的一根
解：在MATLAB命令窗口输入：
>> r=NewtonRoot('3*x^2-exp(x)',3,4)
输出结果：
X=
流程图：
２、编程解决以下科学计算问题。
1)
解：这个方程可用下列步骤来解
fclose(fid);
将以上代码保存为在MATLAB中运行，运行结果如下图所示：
将以下代码保存为文件：
fopen('');
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])
运行结果如下图所示：
if(dfa>dfb)%初始值取两端点导数较大者
root=a-fa/dfa;
else
root=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);%求该点的导数值
（1）用eig函数求出矩阵K-λM的特征值L和特征向量U，U和L满足
（2）在原始方程Mx+Kx=0两端各乘以 及在中间乘以对角矩阵 *U，得
*M* *U* + *K* *x=0
作变量置换 ，得
这是一个对角矩阵方程，即可把它分两个方程：