弦截法非线性方程.doc

科学计算的算法及其举例应用和利用MATLAB自带函数实现 科学计算法：弦截法非线性方程 实例: 采用弦截法求方程lgx+＝2在区间[1，4]上的一个根。 1.算法说明 过两点(a,f(a)),(b,f(b))作一直线，它与x轴有一个交点，记为x1.. 如果f(a)f(x1)0,过两点(a,f(a)),(x1,f(x1))作一直线，它与x轴的交点记为 x2，否则过两点 (b,f(b)),(x1,f(x1))作一直线，它与x轴的交点记为x2 ； 如此下去，直到|xn-xn-1|，就可以认为xn为f(x)=0在区间[a,b]上的一个根。 Xk的递推公式为：且 在MATLAB中编程实现的弦截法的函数为：Secant. 功能：用弦截法求函数在某个区间的一个零点。 调用格式：root=Secant(f,a,b,eps). 其中，f为函数名； a为区间左端点； b为区间左端点； eps为根的精度； root为求出的函数零点。 2．流程图 弦截法流程图 实例流程图： 3．源程序代码 function root=Secant(f,a,b,eps) 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*f20) 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); else root=a-(r1-a)*fa/(fx-fa); end end tol=abs(root-r1); end end 实例源代码： r=Secant(sqrt(x)+log(x)-2,1,4) 运行结果： r = 1.8773 分析单自由度阻尼系统的阻尼系数对其固有振动模态的影响 流程图 2.源程序代码： 3.运行结果： 三．Appollo卫星的运行轨迹 1.流程图 2.源程序代码： 函数文件： 命令文件： 3.运行结果 f(a)f(x1)0?？ YES (a,f(a)),(x1,f(x1))作一直线，它与x轴的交点记为 x2 NO (b,f(b)),(x1,f(x1))作一直线，它与x轴的交点记为x2 YES YES xn为f(x)=0在区间[a,b]上的一个根 输出xn为f(x)=0在区间[a,b]上的一个根 过(a,f(a)),(b,f(b))作一直线，它与x轴有一个交点，记为x1. 输入a,b, NO NO 输出结果 r =1.8773 结束 输入条件：sqrt(x)+log(x)-2 积分区间：（1,4） 开始 结束 绘制x,t二维图象 subplot(1,2,1) j≤10? plot(t,x(j,:)) j=j+2 开始 输入wn， x0， v0 ，tf j=1 zeta=0.1*j wd=wn A=/wd a=atan2(wd*x0,v0+zeta*wn*x0) j≤10? t=0 j=j+1 t≤tf ? t=t+tf/1000 x(j,:)=A** sin(wd*t+a) subplot(1,2,2