抛物线法非线性方程求解

《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 进行插值得到一条抛物线，它有两个根：

,24223C

AC

B B x x -±-+

= 其中

。

)](,,[],[],

,,[),(12012120102x x x x x f x x f B x x x f C x f A -+===

两个根中只取靠近2x 的那个根，即±号取于B 同号， 即

AC

B B B A

x x 4)sgn(22

23-+-

=

（3）用321,,x x x 代替210,,x x x ，重复以上步骤，并有以下递推公式： n

n n

n n n

n 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

%求出的函数零点：root

if(nargin==4)

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;

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,y]=ode45('fun',[1,2,3],1)

x =

1

2

3

y =

1.0000

1.3679

2.1353

>> [x,y]=ode23('fun',[1,2,3],1)

x =

1

2

3

y =

1.0000

1.3677

2.1352

>> [x,y]=ode113('fun',[1,2,3],1) x =

1

2

3

y =

1.0000

1.3679

2.1353