中点公式法和五点公式法求数值微分.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,'-'）
升级会员 