弦截法非线性方程求解.doc

MATLAB程序设计实践课程考核一、,编程实现以下科学计算算法，并举一例应用之。（参考书籍精通MATLAB科学计算，王正林等著，电子工业出版社，2009年）“弦截法非线性方程求解”1、算法说明:弦截法的算法过程如下：（1） 过两点(a,f(a),(b,f(b)作一直线，它与x轴有一个交点，记为x1.（2） 如果f(a)f(x1)0,过两点(a,f(a),(x1,f(x1)作一直线，它与x轴的交点记为 x2，否则过两点 (b,f(b),(x1,f(x1)作一直线，它与x轴的交点记为x2 ；（3） 如此下去，直到|xn-xn-1|，就可以认为xn为f(x)=0在区间a,b上的一个根。（4） Xk的递推公式为：且在MATLAB中编程实现的弦截法的函数为：Secant.功能：用弦截法求函数在某个区间的一个零点。调用格式：root=Secant(f,a,b,eps).其中，f为函数名；a为区间左端点；b为区间右端点；eps为根的精度；root为求出的函数零点。 2、流程图:f(a)f(x1)0) disp(两端点函数值大于0!); return;else tol=1;fa=subs(sym(f),findsym(sym(f),a);fb=subs(sym(f),findsym(sym(f),b);root=a-(b-a)*fa./(fb-fa); %迭代初始值while (toleps) r1=root; fx=subs(sym(f),findsym(sym(f),r1); s=fx*fa; if(s=0) root=r1; else if(s0) root=b-(r1-b)*fb/(fx-fb); %用递推公式2 else root=a-(r1-a)*fa/(fx-fa); %用递推公式1 end end tol=abs(root-r1)endend4、程序应用举例: 采用弦截法求方程lgx+2在区间1，4上的一个根。 解：在MATLAB命令窗口中输入： r=Secant(sqrt(x)+log(x)-2,1,4)输出计算结果为：r = 1.8773由计算结果知方程在lgx+2在区间1，4上的一个根为1.87735、运行流程图:调用Secant进行运算输出结果 r =1.8773结束输入条件：sqrt(x)+log(x)-2积分区间：（1,4）开始6、运行结果:二.分析单自由度阻尼系统的阻尼系数对其固有振动模态的影响。1、算法说明:根据题目意思可理解为在一定范围内取不同的阻尼系数值，根据单自由度阻尼系统的动力学方程，则可以分别算出对应的振动曲线，即振动规律曲线。用matlab进行编程计算，通过调用subplot函数分别绘制x-t二维和三维图形，实现其振动规律曲线，并通过曲线对阻尼振动规律进行分析。2、流程图是否否否结束绘制x,t二维图象subplot(1,2,1)j=1j10?plot(t,x(j,:)j=j+2开始输入wn， x0， v0 ，tfj=1zeta=0.1*jwd=wn1-zeta2A=(wn*x0*zeta+v0)2+(x0*wd)2)/wda=atan2(wd*x0,v0+zeta*wn*x0)j10?t=0j=j+1ttf ?t=t+tf/1000x(j,:)=A*e-zata*wn*t* sin(wd*t+a)subplot(1,2,2),mesh(x)是是3、M文件（U.m）：wn=10; %公共参数 x0=1; %初始位置 v0=0; %初始速度 tf=2; %终点时间 for j=1:10 zeta(j)=0.1*j; %阻尼系数 wd(j)=wn*sqrt(1-zeta(j)2); A=sqrt(wn*x0*zeta(j)+v0)2+(x0*wd(j)2)/wd(j); %振幅A a=atan2(wd(j)*x0,v0+zeta(j)*wn*x0); %相位角 t=0:tf/1000:tf x(j,:)=A*exp(-zeta(j)*wn*t).*sin(wd(j)*t+a); end subplot(1,2,1) %绘制x,t二维图象 for j=1:2:10 plot(t,x(j,:) text(0.3,x(j,151),zeta=,num2str(zeta(j) hold on end xlabel(t),ylabel(x) grid subplot(1,2,2),mesh(x) %绘制对应的三维图形4、在MATLAB命令窗口中运行： U5、运行结果图三. 实验4 已知Appolo卫星的运动轨迹（x,y）满足下面方程：，其中， ，试在初值x(0)=1.2,x(0)=0, y(0)=0, y(0)=-1.04935371下求解，并绘制Appolo卫星轨迹图。1、 算法说明:根据题目意思理解，首先建立目标函数appollo(t,x)给常量muw,lamda及变量r1,r2赋值，令x=x ；x； y； y，则dx=x； x； y ；y；再调用常微分方程函数ode45求出题目给出的微分方程组的数值解；最后调用绘图函数plot绘出x和y的图形，即阿波罗卫星的轨迹图。2、流程图：设置积分的相对误差1e-8调用常微分方程函数ode45求出数值解调用绘图函数plot绘出x对y的图形建立目标函数appolo(t,x)给常量muw,lamda及变量r1,r2赋值令x=x ；x； y； y则dx=x； x； y ；y设定初值x(0)=1.2，x(0)=0y(0)=0,y(0)=-1.04935371积分限定为0，20开始结束3、M文件(函数文件appollo.m)：function dx=appollo(t,x)mu=1/82.45;lamda=1-mu;r1=sqrt(x(1)+mu)2+x(3)2);r2=sqrt(x(1)+lamda)2+x(3)2);dx=x(2);2*x(4)+x(1)-lamda*(x(1)+mu)/r13-mu*(x(1)-lamda)/r23;x(4);-2*x(2)+x(3)-lamda*x(3)/r13-mu*x(3)/r23; % 令x=x ；x； y； y 则dx=x； x； y ；y 4、在MATLAB命令窗口中输入： x0=1.2;0;0;-1.04935371; %x0(i)对应与xi的初值options=odeset(reltol,1e-8); %设置积分的相对误差1e-8tict,y=ode45(appollo,0,20,x0,options); %调用常微分方程函数ode45求出数值解并积分限定为0，20%t是时间点，y的第i列对应xi的值，t和