1、MATLAB 程序设计实践1、编程实现以下科学计算算法，并举一例应用之。“中点公式法和五点公式法求数值微分”解：中点公式法和五点公式法求数值微分如下：例 5-4 ：中点公式法求导数应用实例。采用中点公式法求函数 f= 在 x=4 处的导数。x解：在 MATLAB 命令窗口中输入：df=MidPoint(sqrt(x),4)输出结果为：df=0.2500采用中点公式法求函数 f= 在 x=4 处的导数为 0.25，而导数的精确值也是 0.25x.详见以下：中点公式法流程图：开始If nargin=2,h=0.1 if nargin=3,h=0源代码：function df=MidPoint(fu
2、nc,x0,h)if nargin = 2 判断输入参数个数用中点公式求数值微分结束 输出：h 不能为 0h = 0.1;else if (nargin = 3 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),
3、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,从上面的结果来看，五点公式的精度是很高的。详见以下：五点公式法流程图：If nargin=3,h=0.1 if nargin =4,h=0开始判断输入参数个数用公式一求导
4、用公式二求导用公式三求导用公式四求导用公式五求导结束输出: h 不能为 0 源代码：function df=FivePoint(func,x0,type,h)%,funcx0 %func%x0%type1,2,3,4,5,%h%dfif nargin =3h=0.1; else if (nargin =4return;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
5、 = 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);switch typecase 1,df=(-2
6、5*y0+48*y1-36*y2+16*y3-3*y4)/(12*h);%case 2,df=(-3*y_1-10*y0+18*y1-6*y2+y3)/(12*h);%case 3,df=(y_2-8*y_1+8*y1-y2)/(12*h);%case 4,df=(3*y1+10*y0-18*y_1+6*y_2-y_3)/(12*h);%case 5,df=(25*y0-48*y_1+36*y_2-16*y_3+3*y_4)/(12*h);%end运行结果如下：2、编程解决以下科学计算和工程实际问题。已知阿波罗（Apollo）卫星的运动轨迹 (x,y)满足下列微分方程rxxy32*31*. )
7、(2y321*.其中 = ， =1- ， 试在初45.821*2)(xr2*)(yxr值 x(0)=1.2, , 下进行数值求解，并绘制出阿波罗卫星位置0)(x,049357.)(.y(x,y)的轨迹。解：根据题目选用 MATLAB 代码如下：function dy=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);d
8、y(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,0 2.00,1.2 0 0 -1.04935751)T,Y=ode45(weifen,0 1.26,1.2 0 0 -1.04935751)运行结果：阿波罗卫星位置(x,y)的轨迹图如下：实验图所示是一个跷跷板，两板价较为，左边板长为 1.5m，上面的小孩重 150N,右边板长为 2m,小孩重为 400N.求当跷跷板平衡时，左边木板与水平方向的夹角 的大小。要求先求解析解，然后给出两种解决方案。解：根据力矩平衡求解析解 由图示可有下列关系式：500 1.5 =2 400cos)31cos(解该式得： 738arctni4035即： d46.0两种方法的求解：方案一：采用两步迭代法求解方程：500 1.5 =2 400cos)31cos(两步迭代法的 MATLAB 的代码如下：
