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

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

《中点公式法和五点公式法求数值微分.docx》由会员分享，可在线阅读，更多相关《中点公式法和五点公式法求数值微分.docx（16页珍藏版）》请在冰豆网上搜索。

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

《MATLAB程序设计实践》

1、编程实现以下科学计算算法，并举一例应用之。

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

解：

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

例5-4：

中点公式法求导数应用实例。

采用中点公式法求函数f=

在x=4处的导数。

解：

在MATLAB命令窗口中输入：

>>df=MidPoint（'sqrt（x）',4）

输出结果为：

df=

0.2500

采用中点公式法求函数f=

在x=4处的导数为0.25，而导数的精确值也是0.25

.

详见以下：

中点公式法流程图：

Ifnargin=2,h=0.1ifnargin=3,h=0

源代码：

functiondf=MidPoint（func,x0,h）

ifnargin==2

h=0.1;

elseif（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,从上面的结果来看，五点公式的精度是很高的。

详见以下：

五点公式法流程图：

Ifnargin=3,h=0.1ifnargin=4,h=0

源代码：

functiondf=FivePoint（func,x0,type,h）

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

%º¯ÊýÃû£ºfunc

%Ç󵼵㣺x0

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

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

%µ¼ÊýÖµ£ºdf

ifnargin==3

h=0.1;

elseif（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）;

switchtype

case1,

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）满足下列微分方程

其中

=

，

=1-

，

试在初值x（0）=1.2,

下进行数值求解，并绘制出阿波罗卫星位置（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

=2

400

解该式得：

即：

两种方法的求解：

方案一：

采用两步迭代法求解方程：

500

1.5

=2

400

两步迭代法的MATLAB的代码如下：

源代码：

functionroot=TwoStep（f,a,b,type,eps）

if（nargin==4）

eps=1.0e-4;%eps±íʾµü´ú¾«¶È

end

f1=subs（sym（f）,findsym（sym（f））,a）;

f2=subs（sym（f）,findsym（sym（f））,b）;

if（f1==0）

root=a;

end

if（f2==0）

root=b;

end

if（f1*f2>0）

disp（'两端点函数值乘积大于0!

'）;

return;

else

tol=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（dfa>dfb）

root=a;

else

root=b;

end

while（tol>eps）

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）;

else

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*（r2-r1）/（2*fx2-fx1）;

tol=abs（root-r1）;

end

end

end

解：

在MATLAB命令窗口中输入

r=TwoStep（'500*1.5*cos（x）-2*400*cos（1/3*pi-x）',0,1/3*pi,1）

运行结果如下：

两个小孩所产生力矩随

变化的曲线如下图片：

运用Datacursor工具，得到交点处对应的X值为:

0.4678

也即：

>>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,'-'）