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

《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

.

详见以下：

中点公式法流程图：

If nargin=2,h=0.1 if nargin=3,h=0

源代码：

function df=MidPoint(func,x0,h) if nargin == 2

判断输入参

数个数

用中点公式求数值微分 结束

输出：h 不能为0

h = 0.1;

else if (nargin == 3 && h == 0.0)

disp('h²»ÄÜΪ0£¡');

return;

end

end

y1 = 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.4161

df2=

-0.4161

df3=

-0.4161

df4=

-0.4161

df5=

-0.4161

而函数在f=sin(x)在x=2的导数为cos(2)=-0.4161,从上面的结果来看，五点公式的精度是很高的。

详见以下：

五点公式法流程图：

源代码：

function df=FivePoint(func,x0,type,h)

%Îåµã¹«Ê½·¨,ÇóÈ¡º¯ÊýfuncÔÚx0´¦µ¼Êý ×£ÍòÊÂÈçÒâ

%º¯ÊýÃû£ºfunc

%Ç󵼵㣺x0

%¹«Ê½µÄÐÎʽ£ºtype£¨È¡1,2,3,4,5,£©

%ÀëÉ¢²½³¤£ºh

%µ¼ÊýÖµ£ºdf

if nargin ==3

h=0.1;

else if (nargin ==4&&h==0.0)

disp('h²»ÄÜΪ0');

return;

end

end

y0 = 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);

switch type

case 1,

df=(-25*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)满足下列微分方程

()r

r

x x x y

x 32

*31

*..

)

(2μμμμ--

+-+=

r

r

y

y

y x y 32

31

*.

..

2μμ-

-+-=

其中μ=45

.821

，*

μ=1-μ

221

)(y x r

++=μ ，22*2)(y x r ++=μ 试在初值

x(0)=1.2, 0)0(.

=x , ,04935751.1)0(.

-=y 下进行数值求解，并绘制出阿波罗卫星位置(x,y)