用拋物線與橫軸的交點(diǎn)來逼近函數(shù)f(x)的根

《MATLAB程序設(shè)計(jì)實(shí)踐》課程考核１、編程實(shí)現(xiàn)以下科學(xué)計(jì)算算法，并舉一例應(yīng)用之。（參考書籍《精通ＭＡＬＡＢ科學(xué)計(jì)算》，王正林等著，電子工業(yè)出版社，２００９年）“拋物線法非線性方程求解”拋物線法弦截法其實(shí)是用不斷縮短的直線來近似函數(shù)f（x）的，而拋物線法則是采用三個(gè)點(diǎn)來近似函數(shù)f（x），并且用拋物線與橫軸的交點(diǎn)來逼近函數(shù)f（x）的根。它的算法過程如下：選定初始值x0，x1，x2，并計(jì)算f（x0），f（x1），f（x2）和以下差分：f[x2.x1]=ff[x1,x0]=ff[x2,x1,x0]=f一般取x0=a，x1=b，a<x2<b。注意不要使三點(diǎn)共線。用牛頓插值法對三點(diǎn)（x0，f（x0）），（x1，f（x1））。（x2，f（x2））進(jìn)行插值得一條拋物線，它有兩個(gè)根：x3=x2+-B±其中A=f(x2),C=f[x2,x1,x0]，B=f[x2,x1]+f[x2,x1,x0](x2-x1).兩個(gè)根中只取靠近x2的那個(gè)根，即號取與Ｂ同號，即x3=x2--2A（３）用x1，x2，x3代替ｘ０，ｘ１，ｘ２，重復(fù)以上步驟，并有以下遞推公式：xn+1=xn-2其中AnBn（４）進(jìn)行精度控制。在MATLAB中編程實(shí)現(xiàn)的拋物線法的函數(shù)為：Parabola。功能：用拋物線法求函數(shù)在某個(gè)區(qū)間上的一個(gè)零點(diǎn)。調(diào)用格式：root=Parabola（f，a，b，x，eps）。其中，f為函數(shù)名：a為區(qū)間左端點(diǎn)：b為區(qū)間右端點(diǎn)：x為初始迭代點(diǎn)：eps為根的精度：root為求出的函數(shù)零點(diǎn)。拋物線法的MATLAB程序代碼如下：functionroot=Parabola（f，a，b，x，eps）%拋物線法求函數(shù)f在區(qū)間【a，b】上的一個(gè)零點(diǎn)%函數(shù)名：f%區(qū)間左端點(diǎn)：a%區(qū)間右端點(diǎn)：b%初始迭代點(diǎn)：x%根的精度：eps%求出的函數(shù)零點(diǎn)：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(‘兩端點(diǎn)函數(shù)值乘積大于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個(gè)點(diǎn)t(2)=t(3);t(3)=root;f1=subs(sym(f),findsym(sym(f)),t(1));%計(jì)算3個(gè)點(diǎn)的函數(shù)值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));%計(jì)算3個(gè)差分d2=(f3-f2)/(t(3)-t(2));d3=(d2-d1)/(t(3)-t(1));B=d2+d3*(t(3)-t(2));%計(jì)算算法中的Broot=t(3)-2*f3/(B+sign(B)*sqrt(B^2-4*f3*d3));tol=abs(root-t(3));endend例子：拋物線法求解非線性方程求方程lgx+x=2在區(qū)間【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);進(jìn)行牛頓插值法root=x-2*fx/(B+sign(B)*sqrt(B^2-4*fx*d3));輸入?yún)?shù)a、b、x的值，輸出r=a輸出r=b輸出兩端點(diǎn)函數(shù)值乘積大于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;計(jì)算函數(shù)f（a）、f（b），并令f1=f（a）、f2=f（b）計(jì)算函數(shù)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計(jì)算算法中的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否結(jié)束開始輸入函數(shù)及變量開始輸入函數(shù)及變量r=Parabola(‘sqrt(x)+log(x)-2,1,4,2)調(diào)用Parabola進(jìn)行運(yùn)算輸出結(jié)果r=1.8773結(jié)束實(shí)驗(yàn)2，編程解決以下科學(xué)計(jì)算問題。用ode45，ode23和ode113解下列微分方程：y=x-y，y(0)=1,0<x<3(要求輸出x=1,2,3點(diǎn)的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)得到結(jié)果：輸入（2）：[x,y]=ode23('fun',[1,2,3],1)得到結(jié)果：輸入（3）：[x,y]=ode113('fun',[1,2,3],1)得到結(jié)果：2.M文件：functionf=f(x,y)f=[23;2-1]*y;輸入（1）：ode45(‘fun2’,[0,10],[-2.7,2.8])得到結(jié)果：輸入（2）：ode23('fun2',[0,10],[-2.7,2.8])得到結(jié)果：輸入（3）：ode113('fun2',[0,10],[-2.7,2.8])得到結(jié)果：3.將多階微分方程變?yōu)橐浑A微分方程因此原方程可改寫成：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))得到結(jié)果：輸入（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))輸出結(jié)果：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))輸出結(jié)果：輸入（2）：[t,x]=ode23('fun4',[0,2],[2,1]);plot(t,x(:,1))輸出結(jié)果：輸入（3）：[t,x]=ode113('fun4',[0,2],[2,1]);plot(t,x(:,1))輸出結(jié)果：5.當(dāng)=1時(shí)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)輸出結(jié)果：輸入（2）：[x,y]=ode23('Vanderpol',[0,20],[2,0]);plot(x,y)輸出結(jié)果：輸入（3）：[x,y]=ode113('Vanderpol',[0,20],[2,0]);plot(x,y)輸出結(jié)果：當(dāng)=2時(shí)，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)輸出結(jié)果：輸入（2）：[x,y]=ode23('Vanderpol2',[0,20],[2,0]);plot(x,y)輸出結(jié)果：輸入（3）：[x,y]=ode113('Vanderpol2',[0,20],[2,0]);plot(x,y)輸出結(jié)果：實(shí)驗(yàn)3用ode45和ode15s求解下列剛性方程組，比較計(jì)算速率｛0<x<50M文件functionf=f(x,y)f=[-1000.25999.75;999.75-1000.25]*y+[0.5;0.5];%輸入題目中的參數(shù)在命令窗口輸入：clearode45(‘f’,[050],[1,-1]);ode45clearode15s(‘f’,[050],[1,-1]);ode15s比較計(jì)算效率：M文件t1：tic;[x,y]=ode45('f',[0,50],[1,-1]);t=toc得到結(jié)果：M文件t2：tic;[x,y]=ode15s('f',[0,50],[1,-1]);t=toc得到結(jié)果：即可知ode15s的計(jì)算效率比ode45高很多。實(shí)驗(yàn)4.已知Appolo衛(wèi)星的運(yùn)動軌跡（x,y）滿足下面方程：，其中，，試在初值x(0)=1.2,x’(0)=0,y(0)=0,y’(0)=-1.04935371下求解，并繪制Appolo衛(wèi)星軌跡圖。流程圖設(shè)置積分的相對誤差1e-8設(shè)置積分的相對誤差1e-8調(diào)用常微分方程函數(shù)ode45求出數(shù)值解調(diào)用繪圖函數(shù)plot繪出x對y的圖形建立目標(biāo)函數(shù)appolo(t,x)給常量mu,lamda及變量r1,r2賦值令x=[x；x’；y；y’]則dx=[x’；x’’；y’；y’’]設(shè)定初值x(0)=1.2，x’(0)=0y(0)=0,y’(0)=-1.04935371積分限定為[0，20]開始結(jié)束（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);