中点公式法和五点公式法求数值微分

《MATLAB程序设计实践》1、编程实现以下科学计算算法，并举一例应用之。“中点公式法和五点公式法求数值微分”解：中点公式法和五点公式法求数值微分如下：例5-4：中点公式法求导数应用实例。采用中点公式法求函数f=x在x=4处的导数。解：在MATLAB命令窗口中输入：df=MidPoint('sqrt(x)',4)输出结果为：df=0.2500采用中点公式法求函数f=x在x=4处的导数为0.25，而导数的精确值也是0.25.详见以下：中点公式法流程图：开始Ifnargin=2,h=0.1ifnargin=3,h=0源代码：functiondf=MidPoint(func,x0,h)ifnargin==2判断输入参数个数用中点公式求数值微分结束输出：h不能为0h=0.1;elseif(nargin==3&&h==0.0)disp('h²»ÄÜΪ0£¡');return;endendy1=subs(sym(func),findsym(sym(func)),x0+h);y2=subs(sym(func),findsym(sym(func)),x0-h);df=(y1-y2)/(2*h);运行结果如下：例5-5：五点公式法求导数应用实例。采用五点公式法求函数f=sin(x)在x=2处的导数。解：在MATLAB命令窗口中输入：df1=FivePoint('sin(x)',2,1);df2=FivePoint('sin(x)',2,2);df3=FivePoint('sin(x)',2,3);df4=FivePoint('sin(x)',2,4);df5=FivePoint('sin(x)',2,5);用五种方法得到的结果为：df1=-0.4161df2=-0.4161df3=-0.4161df4=-0.4161df5=-0.4161而函数在f=sin(x)在x=2的导数为cos(2)=-0.4161,从上面的结果来看，五点公式的精度是很高的。详见以下：五点公式法流程图：Ifnargin=3,h=0.1ifnargin=4,h=0开始判断输入参数个数用公式一求导用公式二求导用公式三求导用公式四求导用公式五求导结束输出:h不能为0源代码：functiondf=FivePoint(func,x0,type,h)%Îåµã¹«Ê½·¨,ÇóÈ¡º¯ÊýfuncÔÚx0´¦µ¼Êý×£ÍòÊÂÈçÒâ%º¯ÊýÃû£ºfunc%Ç󵼵㣺x0%¹«Ê½µÄÐÎʽ£ºtype£¨È¡1,2,3,4,5,£©%ÀëÉ¢²½³¤£ºh%µ¼ÊýÖµ£ºdfifnargin==3h=0.1;elseif(nargin==4&&h==0.0)disp('h²»ÄÜΪ0');return;endendy0=subs(sym(func),findsym(sym(func)),x0);y1=subs(sym(func),findsym(sym(func)),x0+h);y2=subs(sym(func),findsym(sym(func)),x0+2*h);y3=subs(sym(func),findsym(sym(func)),x0+3*h);y4=subs(sym(func),findsym(sym(func)),x0+4*h);y_1=subs(sym(func),findsym(sym(func)),x0-h);y_2=subs(sym(func),findsym(sym(func)),x0-2*h);y_3=subs(sym(func),findsym(sym(func)),x0-3*h);y_4=subs(sym(func),findsym(sym(func)),x0-4*h);switchtypecase1,df=(-25*y0+48*y1-36*y2+16*y3-3*y4)/(12*h);%ÓõÚÒ»¸ö¹«Ê½Çóµ¼Êýcase2,df=(-3*y_1-10*y0+18*y1-6*y2+y3)/(12*h);%Óõڶþ¸ö¹«Ê½Çóµ¼Êýcase3,df=(y_2-8*y_1+8*y1-y2)/(12*h);%ÓõÚÈý¸ö¹«Ê½Çóµ¼Êýcase4,df=(3*y1+10*y0-18*y_1+6*y_2-y_3)/(12*h);%ÓõÚËĸö¹«Ê½Çóµ¼Êýcase5,df=(25*y0-48*y_1+36*y_2-16*y_3+3*y_4)/(12*h);%ÓõÚÎå¸ö¹«Ê½Çóµ¼Êýend运行结果如下：2、编程解决以下科学计算和工程实际问题。①已知阿波罗（Apollo）卫星的运动轨迹(x,y)满足下列微分方程rrxxxyx32*31*..)(2rryyyxy3231*...2其中=45.821，*=1-221)(yxr，22*2)(yxr试在初值x(0)=1.2,0)0(.x,,04935751.1)0(.y下进行数值求解，并绘制出阿波罗卫星位置(x,y)的轨迹。①解：根据题目选用MATLAB代码如下：functiondy=weifen(t,y)%编程解决阿波罗（Apollo）卫星的运动轨迹求解器属于变步长的一种，采用Runge-Kutta算法万事如意u=1/82.45;b=1-u;dy=zeros(4,1);r=zeros(2,1);r(1)=sqrt((y(1)+u)^2+(y(3))^2);r(2)=sqrt((y(1)+b)^2+(y(3))^2);dy(1)=y(2);dy(2)=2*dy(3)+y(1)-b*(y(1)+u)/(r(1)^3)-u*(y(1)-b)/(r(2)^3);dy(3)=y(4);dy(4)=-2*dy(1)+y(3)-b*y(3)/(r(1)^3)-u*y(3)/(r(2)^3);解：在MATLAB命令窗口中输入ode45('weifen',[02.00],[1.200-1.04935751])[T,Y]=ode45('weifen',[01.26],[1.200-1.04935751])运行结果：阿波罗卫星位置(x,y)的轨迹图如下：②实验图所示是一个跷跷板，两板价较为，左边板长为1.5m，上面的小孩重150N,右边板长为2m,小孩重为400N.求当跷跷板平衡时，左边木板与水平方向的夹角ɑ的大小。要求先求解析解，然后给出两种解决方案。②解：根据力矩平衡求解析解由图示可有下列关系式：5001.5cos=2400)31cos(解该式得：738arctan738tansin3400cos350即：rad4678.0两种方法的求解：方案一：采用两步迭代法求解方程：5001.5cos=2400)31cos(两步迭代法的MATLAB的代码如下：源代码：functionroot=TwoStep(f,a,b,type,eps)if(nargin==4)eps=1.0e-4;%eps±íʾµü´ú¾«¶È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*f20)disp('两端点函数值乘积大于0!');return;elsetol=1;fun=diff(sym(f));fa=subs(sym(f),findsym(sym(f)),a);fb=subs(sym(f),findsym(sym(f)),b);dfa=subs(sym(fun),findsym(sym(fun)),a);dfb=subs(sym(fun),findsym(sym(fun)),b);if(dfadfb)root=a;elseroot=b;endwhile(toleps)if(type==1)r1=root;fx1=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1);r2=r1-fx1/dfx;fx2=subs(sym(f),findsym(sym(f)),r2);root=r2-fx2/dfx;tol=abs(root-r1);elser1=root;fx1=subs(sym(f),findsym(sym(f)),r1);dfx=subs(sym(fun),findsym(sym(fun)),r1);r2=r1-fx1/dfx;fx2=subs(sym(f),findsym(sym(f)),r2);root=r2-fx2*(r2-r1)/(2*fx2-fx1);tol=abs(root-r1);endendend解：在MATLAB命令窗口中输入r=TwoStep('500*1.5*cos(x)-2*400*cos(1/3*pi-x)',0,1/3*pi,1)运行结果如下：两个小孩所产生力矩随变化的曲线如下图片：运用Datacursor工具，得到交点处对应的X值为:0.4678也即：4678.0x=0:0.0001:1/3*pi;ML=500*1.5*cos(x);MR=400*2*cos(1/3*pi-x);plot(x,ML,'-',x,MR,'-')