用抛物线与横轴的交点来逼近函数f(x)的根

《MATLAB程序设计实践》课程考核１、编程实现以下科学计算算法，并举一例应用之。（参考书籍《精通ＭＡＬＡＢ科学计算》，王正林等著，电子工业出版社，２００９年）“抛物线法非线性方程求解”抛物线法弦截法其实是用不断缩短的直线来近似函数f（x）的，而抛物线法则是采用三个点来近似函数f（x），并且用抛物线与横轴的交点来逼近函数f（x）的根。它的算法过程如下：选定初始值x0，x1，x2，并计算f（x0），f（x1），f（x2）和以下差分：f[x2.x1]=ff[x1,x0]=ff[x2,x1,x0]=f一般取x0=a，x1=b，a<x2<b。注意不要使三点共线。用牛顿插值法对三点（x0，f（x0）），（x1，f（x1））。（x2，f（x2））进行插值得一条抛物线，它有两个根：x3=x2+-B±其中A=f(x2),C=f[x2,x1,x0]，B=f[x2,x1]+f[x2,x1,x0](x2-x1).两个根中只取靠近x2的那个根，即号取与Ｂ同号，即x3=x2--2A（３）用x1，x2，x3代替ｘ０，ｘ１，ｘ２，重复以上步骤，并有以下递推公式：xn+1=xn-2其中AnBn（４）进行精度控制。在MATLAB中编程实现的抛物线法的函数为：Parabola。功能：用抛物线法求函数在某个区间上的一个零点。调用格式：root=Parabola（f，a，b，x，eps）。其中，f为函数名：a为区间左端点：b为区间右端点：x为初始迭代点：eps为根的精度：root为求出的函数零点。抛物线法的MATLAB程序代码如下：functionroot=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));%计算算法中的Broot=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3));tol=abs(root-t(3));endend例子：抛物线法求解非线性方程求方程lgx+x=2在区间【1,4】上的一根。流程图开始否开始否是否是否是否是是f1=0f2=0f1*f2〉0让tol=1；fa=f1，fb=f2；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));输入参数a、b、x的值，输出r=a输出r=b输出两端点函数值乘积大于0!建立3阶0矩阵t，并令t（1）=a，t（2）=b，t（3）=xtol>eps令t(1)=t(2);t(2)=t(3);t(3)=root;计算函数f（a）、f（b），并令f1=f（a）、f2=f（b）计算函数f[t（1）]、f[t（2）]、f[（t（3）]的值，并令d1=(f2-f1)/(t(2)-t(1));d2=(f3-f2)/(t(3)-t(2));d3=(d2-d1)/(t(3)-t(1));Nargin=4令eps=1.0e-4计算算法中的B:B=d2+d3*(t(3)-t(2));令root=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3));tol=abs(root-t(3));输出root否结束开始输入函数及变量开始输入函数及变量r=Parabola(‘sqrt(x)+log(x)-2,1,4,2)调用Parabola进行运算输出结果r=1.8773结束实验2，编程解决以下科学计算问题。用ode45，ode23和ode113解下列微分方程：y=x-y，y(0)=1,0<x<3(要求输出x=1,2,3点的y值)；x=2x+3y,y=2x-y,x(0)=-2.7,y(0)=2.8,0<t<10,作相平面图；y-0.01(y)+2y=sin(t),y(0)=0,y=1,0<t<5,作y的图；2x(t)-5x(t)-3x(t)=45e,x(0)=2,x(0)=1,0<t<2,作x的图；Vanderpol方程y-(y-1)y-y=0,y(0)=2,y(0)=0,0<x<20,=1和2，作相平面图。解：1.M文件：functionf=f(x,y)f=x-y;输入（1）：[x,y]=ode45('fun',[1,2,3],1)得到结果：输入（2）：[x,y]=ode23('fun',[1,2,3],1)得到结果：输入（3）：[x,y]=ode113('fun',[1,2,3],1)得到结果：2.M文件：functionf=f(x,y)f=[23;2-1]*y;输入（1）：ode45(‘fun2’,[0,10],[-2.7,2.8])得到结果：输入（2）：ode23('fun2',[0,10],[-2.7,2.8])得到结果：输入（3）：ode113('fun2',[0,10],[-2.7,2.8])得到结果：3.将多阶微分方程变为一阶微分方程因此原方程可改写成：M文件：functionxp=fun3(t,y)xp=zeros(2,1)xp(1)=y(2)xp(2)=0.01*(y(2))^2-2*y(1)+sin(t);输入（1）：[t,y]=ode45('fun3',[0,5],[0,1]);plot(t,y(:,1))得到结果：输入（2）：[t,y]=ode23('fun3',[0,5],[0,1]);plot(t,y(:,1))输出：输入（3）：[t,y]=ode113('fun3',[0,5],[0,1]);plot(t,y(:,1))输出结果：4.M文件：functionxp=fun4(t,x)xp=zeros(2,1)xp(1)=x(2)xp(2)=(5*x(2)+3*x(1)+45*exp(2*t))/2;输入（1）；[t,x]=ode45('fun4',[0,2],[2,1]);plot(t,x(:,1))输出结果：输入（2）：[t,x]=ode23('fun4',[0,2],[2,1]);plot(t,x(:,1))输出结果：输入（3）：[t,x]=ode113('fun4',[0,2],[2,1]);plot(t,x(:,1))输出结果：5.当=1时M文件：functionxp=Vanderpol(x,y)xp=zeros(2,1)xp(1)=y(2);xp(2)=(y(1)^2-1)*y(2)-y(1);输入（1）[x,y]=ode45('Vanderpol',[0,20],[2,0]);plot(x,y)输出结果：输入（2）：[x,y]=ode23('Vanderpol',[0,20],[2,0]);plot(x,y)输出结果：输入（3）：[x,y]=ode113('Vanderpol',[0,20],[2,0]);plot(x,y)输出结果：当=2时，M文件：functionxp=Vanderpol(x,y)xp=zeros(2,1)xp(1)=y(2);xp(2)=2*(y(1)^2-1)*y(2)-y(1);输入（1）：[x,y]=ode45('Vanderpol2',[0,20],[2,0]);plot(x,y)输出结果：输入（2）：[x,y]=ode23('Vanderpol2',[0,20],[2,0]);plot(x,y)输出结果：输入（3）：[x,y]=ode113('Vanderpol2',[0,20],[2,0]);plot(x,y)输出结果：实验3用ode45和ode15s求解下列刚性方程组，比较计算速率｛0<x<50M文件functionf=f(x,y)f=[-1000.25999.75;999.75-1000.25]*y+[0.5;0.5];%输入题目中的参数在命令窗口输入：clearode45(‘f’,[050],[1,-1]);ode45clearode15s(‘f’,[050],[1,-1]);ode15s比较计算效率：M文件t1：tic;[x,y]=ode45('f',[0,50],[1,-1]);t=toc得到结果：M文件t2：tic;[x,y]=ode15s('f',[0,50],[1,-1]);t=toc得到结果：即可知ode15s的计算效率比ode45高很多。实验4.已知Appolo卫星的运动轨迹（x,y）满足下面方程：，其中，，试在初值x(0)=1.2,x’(0)=0,y(0)=0,y’(0)=-1.04935371下求解，并绘制Appolo卫星轨迹图。流程图设置积分的相对误差1e-8设置积分的相对误差1e-8调用常微分方程函数ode45求出数值解调用绘图函数plot绘出x对y的图形建立目标函数appolo(t,x)给常量mu,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]开始结束（2）源程序代码M文件functiondx=appolo(t,x)u=1/82.45;l=1-u;r1=sqrt((x(1)+u)^2+x(3)^2);r2=sqrt((x(1)+l)^2+x(3)^2);