《MATLAB 程序设计实践》课程考核
---第37-38页
题１
： 编程实现以下科学计算算法，并举一例应用之。（参考书籍《精
通ＭＡTＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９
年）
“牛顿法非线性方程求解”
弦截法本质是一种割线法，它从两端向中间逐渐逼近方程的根；牛顿法本质上是一种切线法，
它从一端向一个方向逼近方程的根，其递推公式为：
nnxx
1
)()('nnxf
xf
初始值可以取)('af和)('bf的较大者，这样可以加快收敛速度。
和牛顿法有关的还有简化牛顿法和牛顿下山法。
在MATLAB中编程实现的牛顿法的函数为：NewtonRoot。
功能：用牛顿法求函数在某个区间上的一个零点。
调用格式：root=NewtonRoot)(```epsbaf
其中，f为函数名；
a
为区间左端点；
b
为区间右端点
eps为根的精度；
root为求出的函数零点。
，
牛顿法的matlab程序代码如下：
function root=NewtonRoot(f,a,b,eps)
%牛顿法求函数f在区间[a,b]上的一个零点
%函数名:f
%区间左端点:a
N
Y
输入参数值
nargin==3
f1==0
eps=1.0e-4
结束
f2==0
f1*f2>0
两端点函数值乘
积大于0 ！
Y
N
dfa>dfb
Y
N
N
N
Y
Y
开始
tol>eps
Y
输出结果
N
流程图：
root=a
root=b
root=a-fa/dfa
root=b-fb/dfb
root=r1-fx/dfx
%区间右端点:b
%根的精度:eps
%求出的函数零点:root
if(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;
end
if (f2==0)
root=b;
end
if (f1*f2>0)
disp('两端点函数值乘积大于0 ！');
return;
else
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);
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); %求该点的导数值
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=3.7331
2、编程解决以下科学计算问题
1） 二自由度可解耦系统的振动模态分析，二自由度振动系统如图所示，其一般方程为：
0)12(2)12(2220221)21(221)21(11........xxKxxcxm
xKxKKxcxccxm
可写成矩阵形式：0...KxxCxM
设C=0，即无阻尼情况，则系统可解耦为两种独立的振动模态。
输入参数值
结束
开始
流程图：
构成参数矩阵
设定初始条件
i=1:round(tf/dt)+1
循环计算矩阵指数
绘制图象
系统解耦的振动模态的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)+1
t(i)=dt*(i-1);
y(:,i)=expm(A*t(i))*y0;%循环计算矩阵指数
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命令窗口中输入数据：
即可得出该振动的两种模态
2)用GUI方式解下列PDE：
;0,4sin0),3(30,40,030402222yyxxuxu
uyyu
yx
yux
u
解：第一步，在MATLAB命令窗口输入命令pdetool打开工具箱，调整x坐标范围为[0 4],y
坐标范围为[0 3].通过options选项的Axes Linits设定如下图所示。
第二步，设定矩形区域。点击工具箱栏中的按钮“”，拖动鼠标画出一矩形，并双击该
矩形，设定矩形大小，如下图所示。
第三步，设边界条件。点击工具栏中的按钮“”，并双击矩形区域的相应的边线在弹出
的对话框中设定边界条件。如下图所示，分别为各边框的边界条件。
第四步，设定方程。单机工具栏中的按钮“”，在PDE模式下选择方程类型，如下图
所示，并在其中设定参数。
第五步，单击工具栏中的按钮“”，拆分区域为若干子区域，如下图所示。
第六步，单击工具栏中的按钮“”，将子区域细化，从而保证结果更精确，如下图所
示。
第七步，单击工具栏中的按钮“”，设置所画曲线的特性，如下图所示，并作出其解
的三维图。
第八步，单击图2-8在标出的“plot”按钮，或单击工具栏中的按钮“”，可作出解的
三维图，如下图所示。
简要流程图：
开始
绘制要求区域图
设置边界条件
设置方程参数
剖分网格
作出解的三维图
结束