《MATLAB 程序设计实践》课程考核抛物线法非线性方程求解算法说明：（1）选定初始值210,,x x x ，并计算)(),(),(210x f x f x f 和以下差分： ],[12x x f =1212)()(x x x f x f --10101)()(],[x x x f x f x x f --=20112012],[],[],,[x x x x f x x f x x x f --=一般取b x a b x a x <<==210,,。
注意不要使三点共线。
（2）用牛顿插值法对三点))(,()),(,()),(,(221100x f x x f x x f x 进行插值得到一条抛物线，它有两个根：,24223CACB B x x -±-+= 其中。
)](,,[],[],,,[),(12012120102x x x x x f x x f B x x x f C x f A -+===两个根中只取靠近2x 的那个根，即±号取于B 同号， 即ACB B B Ax x 4)sgn(2223-+-=（3）用321,,x x x 代替210,,x x x ，重复以上步骤，并有以下递推公式： nn nn n nn n C A B B B A x x 4)sgn(221-+-=+，其中。
)](,,[],[],,,[),(121121-------+===n n n n n n n n n n n n n n x x x x x f x x f B x x x f C x f A（4）进行精度控制。
算法流程图：成立抛物线法的MATLAB程序代码如下：function root=Parabola（f，a，b，x，eps）%抛物线法求函数 f在区间【a，b】上的一个零点%函数名：f%区间左端点： a%区间右端点：b%初始迭代点：x%根的精度：eps%求出的函数零点：rootif(nargin==4)eps=1.0e-4;endf1=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;fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),a);fx=subs(sym(f),findsym(sym(f)),x);d1=(fb-fa)/(b-a);d2=(fx-fb)/(x-b);d3=(f2-f1)/(x-a);B=d2+d3*(x-b);root=x-2*fx/(B+sign(B)*sqrt(B^2-4*fx*d3)); t=zeros(3); t(1)=a; t(2)=b; t(3)=x; while(tol>eps)t(1)=t(2); %保存3个点 t(2)=t(3); t(3)=root;f1=subs(sym(f),findsym(sym(f)),t(1)); %计算3个点的函数值 f2=subs(sym(f),findsym(sym(f)),t(2)); f3=subs(sym(f),findsym(sym(f)),t(3));d1=(f2-f1)/(t(2)-t(1)); %计算3个差分 d2=(f3-f2)/(t(3)-t(2)); d3=(d2-d1)/(t(3)-t(1));B=d2+d3*(t(3)-t(2)); %计算算法中的B root=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3)); tol=abs(root-t(3)); end end举例应用：采用抛物线法求方程2lg =+x x 在区间[1，4]上的一个根。
流程图：解：在MA TLAB 命令窗口中输入： >> r=Parabola('sqrt(x)+log(x)-2',1,4,2) r =1.8773实验1用ode45、ode23、ode113解下列微分方程 流程图 编辑相应方程的M 文件是（1）y’=x-y，y(0)=1，0<x<3 (要求输出x=1、2、3点的y值); >> [x,y]=ode45('fun',[1,2,3],1)x =123y =1.00001.36792.1353>> [x,y]=ode23('fun',[1,2,3],1)x =123y =1.00001.36772.1352>> [x,y]=ode113('fun',[1,2,3],1) x =123y =1.00001.36792.1353（2）x’=2x+3y,y’=2x-y,x(0)=-2.7,y(0)=2.8,0<t<10,作相平面图；>>ode45(‘fun2’,[0,10],[-2.7,2.8])>> ode23('fun2',[0,10],[-2.7,2.8])>> ode113('fun2',[0,10],[-2.7,2.8])（3）50,1)0(',0)0(),sin(2)'(01.0''2<<===+-t y y t y y y ,作y 的图； 解：将多阶微分方程变为一阶微分方程⎪⎩⎪⎨⎧==211y dt dy y y因此原方程可改写成：)sin(201.0''122221t y y y y y +-==编辑M 文件fun3.m计算以及绘图程序如下： ode45:>> [t,y]=ode45('fun3',[0,5],[0,1]); >> plot(t,y(:,1))得到y 的图线如下：y 值如下：ode23:>> [t,y]=ode23('fun3',[0,5],[0,1]); >> plot(t,y(:,1))得到y 值的图如下：y值如下：Ode113:>> [t,y]=ode113('fun3',[0,5],[0,1]); >> plot(t,y(:,1))得到y值的图如下：Y值如下：（4）20,1)0(',2)0(,45)(3)('5)(''22<<===--txxetxtxtx t，作x的图；解：ode45：>> [t,x]=ode45('fun4',[0,2],[2,1]);plot(t,x(:,1))Ode23:>> [t,x]=ode23('fun4',[0,2],[2,1]);plot(t,x(:,1))Ode113：>> [t,x]=ode113('fun4',[0,2],[2,1]);plot(t,x(:,1))（5）Vanderpol方程当mu=2时编辑M文件Vanderpol2.mOde45:>> [x,y]=ode45('Vanderpol2',[0,20],[2,0]);plot(x,y)Ode23：>> [x,y]=ode23('Vanderpol2',[0,20],[2,0]);plot(x,y)Ode113:>> [x,y]=ode113('Vanderpol2',[0,20],[2,0]);plot(x,y)当mu=1时编辑M文件Vanderpol.mOde45：>> [x,y]=ode45('Vanderpol',[0,20],[2,0]);plot(x,y)Ode23:>> [x,y]=ode23('Vanderpol',[0,20],[2,0]);plot(x,y)Ode113：>> [x,y]=ode113('Vanderpol',[0,20],[2,0]);plot(x,y)实验3 分别用ode45流程图：解：在命令窗口输入：>>ode45(‘f’,[0,50],[1,-1])>>ode15s(‘f’,[0,50],[1,-1])编制命令文件bj1.m和bj2.m分别运行得>> bj1t =13.2500>> bj2t =0.3900可以看出ode15s的计算效率比ode45高很多。
实验4流程图：解：将二阶微分方程转化成一阶微分方程：⎪⎪⎪⎩⎪⎪⎪⎨⎧====)4()3()2()1(x dt dy x y x dt dx x x源程序代码：function dx=appolo(t,x)mu=1/82.45;lamda=1-mu;x=x(1);y=x(3);r1=sqrt((x(1)+mu)^2+x(3)^2);r2=sqrt((x(1)+lamda)^2+x(3)^2);dx/dt=x(2);dx @/dt 2=2*x(4)+x(1)-lamda*(x(1)+mu)/r1^3-mu*(x(1)-lamda)/r2^3;dy/dt=x(4);dy @/dt 2=-2*x(2)+x(3)-lamda*x(3)/r1^3-mu*x(3)/r2^3];在命令窗口输入如下程序：>>x0=[1.2;0;0;-1.04935371]; （设定边间条件：微分方程的条件，将边界条件存放在一个4行1列的矩阵x0中）>>options=odeset('reltol',1e-8);（options 是积分参数设置条件，误差条件，为ode45设置相对误差：1e-8，reltol 相对误差）>>[t,y]=ode45('appolo',[0,20],x0,options); （ode 是求常微分方程的函数）>>plot(y(:,1),y(:,3)) （plot 是绘制二维图形的函数，用来画卫星轨迹，以y 矩阵第一列为x 轴，以其第三列为y 轴画图）>>title('Appolo 卫星运动轨迹')（图的标题）>>xlabel('X')（X 轴）>>ylabel('Y')（Y 轴）运行结果：。