这个问题可以作为符号运算和数值运算结合的很好的一个例子：利用隐函数求导公式对f
进行一二阶求导，然后利用solve得到B关于t的一、二阶导数的符号表达式，然后再利用eval函数转化成数值表达式：整个代码如下：
1.
2.syms A B ;
3.r=1;beta=pi/5;rho=2;
4.%f的符号表达式
5.f=(r*sin(A)-rho*sin(3*B))*(cos(beta)-sin(beta)*cos(3*B)-...
6. sin(beta)*sin(3*B)*tan(B))-(r*cos(A)-rho*cos(beta)*cos(3*B))*
tan(B);
7.%将A,B分别用90*t和B(t)替换，为的是好利用符号diff函数来求对B关于
t的隐函数F求导
8. F = subs(f,{'A','B'},{'90*t','B(t)'});
9.dFt = diff(F,'t');%一阶导数
10.%将diff(B(t), t)用dBt替换，为的是下一步方便用solve求解diff(B(t), t)
的表达式
11.dFt = subs(dFt,'diff(B(t), t)','dBt');
12.dBt = solve(dFt,'dBt');%得到B关于t的一阶导数的表达式
13.%将dBt用dBt(t)替换，为的是告诉MATLAB,dBt是关于t的函数，能够进一步
求导
14.dFt_ = subs(dFt,'dBt','dBt(t)');
15.ddFt = diff(dFt_,'t');%二阶导数
16.%替换'diff(dBt(t), t)','diff(B(t), t)'，方便求解ddBt的表达式
17.ddFt = subs(ddFt,{'diff(dBt(t), t)','diff(B(t),
t)'},{'ddBt','dBt(t)'});
18.ddBt = solve(ddFt,'ddBt');%求解B关于t的二阶导数的表达式
19.B = @(t) fzero(@(B) (r*sin(90*t)-rho*sin(3*B))*(cos(beta)-...
20. sin(beta)*cos(3*B)-sin(beta)*sin(3*B)*tan(B))-...
21. (r*cos(90*t)-rho*cos(beta)*cos(3*B))*tan(B),1);%B关于t的函数
22.eval(['dBt = @(t) ',char(dBt),';' ])%利用eval函数将符号dBt的表达式
转化为数值函数
23.eval(['ddBt = @(t) ',char(ddBt),';' ])
24.R = 1;
25.C = @(t) R*cos(90*t)/tan(B(t))+sin(90*t);%C的表达式
26.t = 0.2:0.1:2;
27.plot(t,arrayfun(@(T) C(T),t) )%画C关于t的图
28.
复制代码
需要说明的是得到B的函数句柄B(t)后我们可以利用导数的定义来近似表达式dBt，和ddBt，这样的优点是速度快，但是不精确。

上述得到的dBt,ddBt,较为精确，但是计算量比较大。

1.
2.>> (B(1)-B(1.00001))/-0.00001
3.ans =
4.-8.656343*********
5.>> dBt(1)
6.ans =
7.-8.664398751230884
8.>> (dBt(1)-dBt(1.00001))/-0.00001
9.ans =
10. 1.611166395037067e+003
11.>> ddBt(1)
12.ans =
13. 1.610613426031270e+003
14.
复制代码
C关于t的图形。