中点公式法和五点公式法求数值微分Word格式文档下载.docx
《中点公式法和五点公式法求数值微分Word格式文档下载.docx》由会员分享，可在线阅读，更多相关《中点公式法和五点公式法求数值微分Word格式文档下载.docx（16页珍藏版）》请在冰豆网上搜索。
源代码：
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处的导数。
df1=FivePoint（'
sin（x）'
2,1）;
df2=FivePoint（'
2,2）;
df3=FivePoint（'
2,3）;
df4=FivePoint（'
2,4）;
df5=FivePoint（'
2,5）;
用五种方法得到的结果为：
df1=
-0.4161
df2=
df3=
df4=
df5=
而函数在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）
0'
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）;
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（'
[01.26],[1.200-1.04935751]）
运行结果：
阿波罗卫星位置（x,y）的轨迹图如下：
②实验图所示是一个跷跷板，两板价较为，左边板长为1.5m，上面的小孩重150N,右边板长为2m,小孩重为400N.求当跷跷板平衡时，左边木板与水平方向的夹角ɑ的大小。
要求先求解析解，然后给出两种解决方案。
②解：
根据力矩平衡求解析解
由图示可有下列关系式：
500
1.5
=2
400
解该式得：
即：
两种方法的求解：
方案一：
采用两步迭代法求解方程：
500
两步迭代法的MATLAB的代码如下：
functionroot=TwoStep（f,a,b,type,eps）
if（nargin==4）
eps=1.0e-4;
%eps±
í
¾
ü
´
ú
f1=subs（sym（f）,findsym（sym（f））,a）;
f2=subs（sym（f）,findsym（sym（f））,b）;
if（f1==0）
root=a;
if（f2==0）
root=b;
if（f1*f2>
0）
两端点函数值乘积大于0!
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）
else
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）;
root=r2-fx2*（r2-r1）/（2*fx2-fx1）;
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,'
）
升级会员 