隐函数方程的求解[有代码]

taxsh · 发表于 2010-4-16 19:33:27

本帖最后由 messenger 于 2010-4-19 23:23 编辑

(r*sin(A)-ρsin(3B))*(cosβ-sinβcos(3B)-sinβsin(3B)tan(B))=(r*cos(A)-ρcosβcos(3B))*tan(B)。已知：r，ρ，β，求：用A表示B的表达式。即求出B等于A的表达式，（B为未知数）或者用什么软件可以解出来吗？MATLAB我不会求，可以用Matlab做吗？有会做的，请教。谢谢！
syms A B;r=1;beta=pi/5;rho=2;f=(r*sin(A)-rho*sin(3*B))*(cos(beta)-sin(beta)*cos(3*B)-sin(beta)*sin(3*B)*tan(B))-(r*cos(A)-rho*cos(beta)*cos(3*B))*tan(B)
>>B=solve(f,B)，这个不好算，出错。

lin2009 · 发表于 2010-4-17 09:52:18

用A表示B的表达式。即求出B等于A的表达式，（B为未知数）。
意思不够简洁明确，越说越不清楚。实际上就是求变量B的表达式。

Matlab不擅长符号运算，可以用Maple，mathematica等软件求解。
解析式很复杂：
arctan(sin(A)*(-2*sin((1/5)*pi)*ComExpr^2+ComExpr*cos((1/5)*pi)+sin((1/5)*pi))/(-16*sin((1/5)*pi)*ComExpr^4+12*sin((1/5)*pi)*ComExpr^2+4*ComExpr*cos((1/5)*pi)-2*sin((1/5)*pi)+cos(A)), ComExpr)
ComExpr = RootOf(256*_Z^10*sin((1/5)*Pi)^2-640*sin((1/5)*Pi)^2*_Z^8-128*_Z^7*sin((1/5)*Pi)*cos((1/5)*Pi)+(592*sin((1/5)*Pi)^2-32*sin((1/5)*Pi)*cos(A))*_Z^6+224*sin((1/5)*Pi)*_Z^5*cos((1/5)*Pi)+(16+56*sin((1/5)*Pi)*cos(A)-268*sin((1/5)*Pi)^2-4*sin((1/5)*Pi)^2*cos(A)^2)*_Z^4+(8*cos(A)*cos((1/5)*Pi)-116*cos((1/5)*Pi)*sin((1/5)*Pi)+4*cos((1/5)*Pi)*sin((1/5)*Pi)*cos(A)^2)*_Z^3+(-28*sin((1/5)*Pi)*cos(A)-15+63*sin((1/5)*Pi)^2+5*sin((1/5)*Pi)^2*cos(A)^2)*_Z^2+(18*cos((1/5)*Pi)*sin((1/5)*Pi)-2*cos((1/5)*Pi)*sin((1/5)*Pi)*cos(A)^2-8*cos(A)*cos((1/5)*Pi))*_Z+4*sin((1/5)*Pi)*cos(A)-cos(A)^2-sin((1/5)*Pi)^2*cos(A)^2-3*sin((1/5)*Pi)^2)
RootOf(...)表示变量_Z一元十次方方程的根。

rocwoods · 发表于 2010-4-17 11:23:39

这样的问题最好还是用数值计算来实现。隐函数和显式函数对于计算机来说没有本质不同，即都是给定一个自变量，计算机通过计算得到因变量，都是通过计算机中的函数实现的。本例可以利用匿名函数+fzero函数来实现：

r=1;beta=pi/5;rho=2;
B_fun = @(A) fzero(@(B) (r*sin(A)-rho*sin(3*B))*(cos(beta)-sin(beta)*cos(3*B)-sin(beta)*sin(3*B)*tan(B))-(r*cos(A)-rho*cos(beta)*cos(3*B))*tan(B),1);
B_fun(1)
B_fun(2)
B_fun(3)
ans =
1.3102
ans =
1.1369
ans =
1.1431

复制代码

更多MATLAB表示匿名函数的用法请参考下面帖子中第6楼：
http://forum.simwe.com/viewthread.php?tid=791513&highlight=

taxsh · 发表于 2010-4-19 14:59:00

本帖最后由 taxsh 于 2010-4-19 15:04 编辑

2# lin2009

兄弟能告诉我你这个解法是在哪个软件里面编程的吗？能把编的程序给我吗？万分感谢！

taxsh · 发表于 2010-4-19 15:05:06

谢谢兄弟的帮忙，我需要的是表达式

rocwoods · 发表于 2010-4-19 16:01:03

本帖最后由 rocwoods 于 2010-4-19 16:04 编辑

谢谢兄弟的帮忙，我需要的是表达式
taxsh 发表于 2010-4-19 15:05

为什么非得要显式表达式？本问题那么复杂的显式表达式反而不好操作，不是前面已经说了吗：
隐表达式和显式表达式对于计算机来说没有本质不同，即都是给定一个自变量，计算机通过计算得到因变量，都是通过计算机中的函数实现的。你完全可以把上面得到的B_fun当成其表达式！好比f(x) = 1+x,的f(x)一样。

taxsh · 发表于 2010-4-19 16:38:06

6# rocwoods

谢谢你的热心帮助！我实在是没办法解决此问题，所以想请你帮忙解决一下,谢谢！上面的式子中A=90*t，t为时间。想求出：当A为自变量，B为应变量时，B关于时间t的函数以及求一阶二阶导数，还有C=R*cos(A)/tan(B)+sin(A),想画出C关于t的函数图像。

rocwoods · 发表于 2010-4-19 18:37:53

本帖最后由 rocwoods 于 2010-4-19 18:40 编辑

你这个问题可以作为符号运算和数值运算结合的很好的一个例子：利用隐函数求导公式对f进行一二阶求导，然后利用solve得到B关于t的一、二阶导数的符号表达式，然后再利用eval函数转化成数值表达式：整个代码如下：

syms A B ;
r=1;beta=pi/5;rho=2;
%f的符号表达式
f=(r*sin(A)-rho*sin(3*B))*(cos(beta)-sin(beta)*cos(3*B)-...
sin(beta)*sin(3*B)*tan(B))-(r*cos(A)-rho*cos(beta)*cos(3*B))*tan(B);
%将A,B分别用90*t和B(t)替换，为的是好利用符号diff函数来求对B关于t的隐函数F求导
F = subs(f,{'A','B'},{'90*t','B(t)'});
dFt = diff(F,'t');%一阶导数
%将diff(B(t), t)用dBt替换，为的是下一步方便用solve求解diff(B(t), t)的表达式
dFt = subs(dFt,'diff(B(t), t)','dBt');
dBt = solve(dFt,'dBt');%得到B关于t的一阶导数的表达式
%将dBt用dBt(t)替换，为的是告诉MATLAB,dBt是关于t的函数，能够进一步求导
dFt_ = subs(dFt,'dBt','dBt(t)');
ddFt = diff(dFt_,'t');%二阶导数
%替换'diff(dBt(t), t)','diff(B(t), t)'，方便求解ddBt的表达式
ddFt = subs(ddFt,{'diff(dBt(t), t)','diff(B(t), t)'},{'ddBt','dBt(t)'});
ddBt = solve(ddFt,'ddBt');%求解B关于t的二阶导数的表达式
B = @(t) fzero(@(B) (r*sin(90*t)-rho*sin(3*B))*(cos(beta)-...
sin(beta)*cos(3*B)-sin(beta)*sin(3*B)*tan(B))-...
(r*cos(90*t)-rho*cos(beta)*cos(3*B))*tan(B),1);%B关于t的函数
eval(['dBt = @(t) ',char(dBt),';' ])%利用eval函数将符号dBt的表达式转化为数值函数
eval(['ddBt = @(t) ',char(ddBt),';' ])
R = 1;
C = @(t) R*cos(90*t)/tan(B(t))+sin(90*t);%C的表达式
t = 0.2:0.1:2;
plot(t,arrayfun(@(T) C(T),t) )%画C关于t的图

复制代码

需要说明的是得到B的函数句柄B(t)后我们可以利用导数的定义来近似表达式dBt，和ddBt，这样的优点是速度快，但是不精确。上述得到的dBt,ddBt,较为精确，但是计算量比较大。

>> (B(1)-B(1.00001))/-0.00001
ans =
-8.656343840729441
>> dBt(1)
ans =
-8.664398751230884
>> (dBt(1)-dBt(1.00001))/-0.00001
ans =
1.611166395037067e+003
>> ddBt(1)
ans =
1.610613426031270e+003

复制代码

C关于t的图形

taxsh · 发表于 2010-4-20 09:03:06

8# rocwoods

为什么在运行的时候会出现以下错误？我用的是Matlab7.0，是不是版本的问题？谢谢
??? Error using ==> solve
Error, (in fsolve) t is in the equation, and is not solved for
Error in ==> sym.solve at 49
[varargout{1:max(1,nargout)}] = solve(S{:});

taxsh · 发表于 2010-4-20 09:13:41

4# taxsh

到底怎么才能获得邀请码呢？我真的是很苦恼！同学有个但是用起来很不方便，我实在是觉得这么论坛相当的不错！我也很需要这个网站！大家谁能给小弟一个邀请码？谢谢了！

lin2009 · 发表于 2010-4-20 09:52:46

4# taxsh
是采用maple的推导的，程序如下：
r := 1;
beta := pi/5;
rho := 2;
f := (r*sin(A) - rho*sin(3*B))*(cos(beta) - sin(beta)*cos(3*B) - sin(beta)*sin(3*B)*tan(B)) - (r*cos(A) - rho*cos(beta)*cos(3*B))*tan(B)
S := solve(f,B)

再将相同的项提出来，用ComExpr表示。

赞同 6# rocwoods 的观点。从ComExpr的表达可以看出，这是一元十次方程的根，没有显性的表达式（没记错的话，一元三次、四次方程的根有显性的表达式，但也非常复杂），其本身就是一个隐性表达式。

闲云孤鹤 · 发表于 2010-4-20 15:13:07

谢谢朋友们的热心帮助！

闲云孤鹤 · 发表于 2010-4-20 16:12:00

我一步一步的检查，发现这个没替换，您帮我看看“diff(B(t), t)用dBt替换”，没有成功，是不是这个导致dBt = solve(dFt,'dBt')出错，提示这个：??? Error using ==> solve
Error, (in fsolve) t is in the equation, and is not solved for

Error in ==> sym.solve at 49
[varargout{1:max(1,nargout)}] = solve(S{:});
您帮我看看
dFt = diff(F,'t')
dFt =
(90*cos(90*t)-6*cos(3*B(t))*diff(B(t),t))*(910872158600853/1125899906842624-2647149443198255/4503599627370496*cos(3*B(t))-2647149443198255/4503599627370496*sin(3*B(t))*tan(B(t)))+(sin(90*t)-2*sin(3*B(t)))*(7941448329594765/4503599627370496*sin(3*B(t))*diff(B(t),t)-7941448329594765/4503599627370496*cos(3*B(t))*diff(B(t),t)*tan(B(t))-2647149443198255/4503599627370496*sin(3*B(t))*(1+tan(B(t))^2)*diff(B(t),t))-(-90*sin(90*t)+2732616475802559/562949953421312*sin(3*B(t))*diff(B(t),t))*tan(B(t))-(cos(90*t)-910872158600853/562949953421312*cos(3*B(t)))*(1+tan(B(t))^2)*diff(B(t),t)

>> dFt = subs(dFt,'diff(B(t), t)','dBt')%没有将diff(B(t), t)用dBt替换，为什么没有替换呢？

dFt =

(90*cos(90*t)-6*cos(3*B(t))*diff(B(t),t))*(910872158600853/1125899906842624-2647149443198255/4503599627370496*cos(3*B(t))-2647149443198255/4503599627370496*sin(3*B(t))*tan(B(t)))+(sin(90*t)-2*sin(3*B(t)))*(7941448329594765/4503599627370496*sin(3*B(t))*diff(B(t),t)-7941448329594765/4503599627370496*cos(3*B(t))*diff(B(t),t)*tan(B(t))-2647149443198255/4503599627370496*sin(3*B(t))*(1+tan(B(t))^2)*diff(B(t),t))-(-90*sin(90*t)+2732616475802559/562949953421312*sin(3*B(t))*diff(B(t),t))*tan(B(t))-(cos(90*t)-910872158600853/562949953421312*cos(3*B(t)))*(1+tan(B(t))^2)*diff(B(t),t)
dBt = solve(dFt,'dBt')
??? Error using ==> solve
Error, (in fsolve) t is in the equation, and is not solved for

Error in ==> sym.solve at 49
[varargout{1:max(1,nargout)}] = solve(S{:});

rocwoods · 发表于 2010-4-20 17:06:05

换个新的版本试一下？我用的R2009a，没有问题。

wanglu · 发表于 2010-4-20 20:39:42

用Forcal求解的代码：

IMSL::ClearImslErr(); //清空ISML错误输出文件“ImslErr.txt”的内容
AB(B:r,ρ,β,pi:A)= pi=3.141592653589793, r=1, β=pi/5, ρ=2, (r*sin(A)-ρ*sin(3*B))*(cos(β)-sin(β)*cos(3*B)-sin(β)*sin(3*B)*tan(B))-[(r*cos(A)-ρ*cos(β)*cos(3*B))*tan(B)];
B(_A:x,i:A)= //由A求B的函数，将验证结果的准确性
{
A=_A, x=1, i=XSLSF::pqrt[HFor("AB"),&x,1e-6],
if[i==10, printff{"\r\nt={1,r}, A={2,r}, B={3,r}, AB={4,r}, 未满足精度要求，返回值仅作参考\r\n",A/90,A,x,AB(x)}],
x
};
C(t:R,A,B)= R=1, A=90*t, B=B(A), R*cos(A)/tan(B)+sin(A); //定义函数C
Bt(t)= B(90*t); //定义B关于t的函数
Bt1(t)= IMSL::DERIV[HFor("Bt"),1,t,0.001,0.01]; //定义Bt关于t的一阶导数
Bt2(t)= IMSL::DERIV[HFor("Bt"),2,t,0.001,0.01]; //定义Bt关于t的二阶导数
//以下进行计算
B[1];
B[2];
B[3];
(Bt(1)-Bt(1.00001))/[-0.00001];
Bt1[1];
Bt2[1];

复制代码

结果：

1.310178796154323
1.136888463437493
1.143092501544323
-8.656343840729441
-8.66439848538108
1610.613441447826

欢迎访问：http://www.forcal.net/

wanglu · 发表于 2010-4-20 20:47:27

用FcCurve绘制的函数C的图形。
Forcal代码：

AB(B:r,ρ,β,pi:A)= pi=3.141592653589793, r=1, β=pi/5, ρ=2, (r*sin(A)-ρ*sin(3*B))*(cos(β)-sin(β)*cos(3*B)-sin(β)*sin(3*B)*tan(B))-[(r*cos(A)-ρ*cos(β)*cos(3*B))*tan(B)];
B(_A:x,i:A)= //由A求B的函数，将验证结果的准确性
{
A=_A, x=1, i=XSLSF::pqrt[HFor("AB"),&x,1e-6],
if[i==10, printff{"\r\nt={1,r}, A={2,r}, B={3,r}, AB={4,r}, 未满足精度要求，返回值仅作参考\r\n",A/90,A,x,AB(x)}],
x
};
_C(t:R,A,B)= R=1, A=90*t, B=B(A), R*cos(A)/tan(B)+sin(A); //定义函数C

复制代码

闲云孤鹤 · 发表于 2010-4-21 20:46:55

14# rocwoods
我运行的那种错误是版本问题吗？我正在下载一个matlabR2009a好像4.05G.听说安装完要13G，是吗？我电脑是老机器，每个盘只有10G，看来得想办法解决啊！

wanglu · 发表于 2010-4-21 20:55:29

这种含三角函数的方程一般是多解的，解与初值有关，将代码中如下句：
A=_A, x=1, i=XSLSF::pqrt[HFor("AB"),&x,1e-6],

改为：
A=_A, x=0.1, i=XSLSF::pqrt[HFor("AB"),&x,1e-6],

其中x为解方程初值。

得到如下图（1000个绘图点）：

同时输出：

t=0.23599999999999993, A=21.239999999999995, B=-0.10699176696776018, AB=6.044367112320187e-002, 未满足精度要求，返回值仅作参考
t=0.25399999999999995, A=22.859999999999996, B=0.6462678324757678, AB=-1.5476520027029841, 未满足精度要求，返回值仅作参考
t=0.25579999999999997, A=23.021999999999998, B=0.13895573563177144, AB=-0.11853262240858076, 未满足精度要求，返回值仅作参考
t=0.39260000000000178, A=35.33400000000016, B=2.3141538447383634, AB=-3.5848958757671654, 未满足精度要求，返回值仅作参考
t=0.44660000000000249, A=40.194000000000223, B=-0.10821562374783429, AB=3.395683060775545e-002, 未满足精度要求，返回值仅作参考
t=0.46460000000000279, A=41.814000000000249, B=-0.10564977140091936, AB=-0.26969909677125903, 未满足精度要求，返回值仅作参考
t=0.46640000000000276, A=41.976000000000248, B=-1.1362423444262002, AB=8.8817841970012523e-016, 未满足精度要求，返回值仅作参考
t=0.53120000000000356, A=47.80800000000032, B=0.12504694031017638, AB=-3.4111443302274003e-002, 未满足精度要求，返回值仅作参考
t=0.54380000000000372, A=48.942000000000334, B=-0.32845957287961941, AB=-2.718588160144364e-009, 未满足精度要求，返回值仅作参考
t=0.58340000000000425, A=52.506000000000384, B=50.014331368132943, AB=0.1344754683491316, 未满足精度要求，返回值仅作参考
t=0.60320000000000451, A=54.288000000000409, B=0.98271464846136847, AB=-2.8327782756375131, 未满足精度要求，返回值仅作参考
t=0.65180000000000515, A=58.662000000000461, B=-24.690782715089522, AB=-7.8825834748386114e-015, 未满足精度要求，返回值仅作参考
t=0.67340000000000544, A=60.606000000000492, B=-9.9847751317575728e-002, AB=-0.26256920139773843, 未满足精度要求，返回值仅作参考
t=0.72560000000000613, A=65.304000000000556, B=3.0254945449961754, AB=7.9741768743701869e-014, 未满足精度要求，返回值仅作参考
t=0.74000000000000621, A=66.600000000000563, B=0.18584616828066736, AB=-3.6074498433307789e-003, 未满足精度要求，返回值仅作参考
t=0.75620000000000653, A=68.05800000000059, B=-0.27424315553824841, AB=1.2863876630575533e-012, 未满足精度要求，返回值仅作参考
t=0.86420000000000796, A=77.778000000000716, B=-0.27595415917319976, AB=0.10993289340065593, 未满足精度要求，返回值仅作参考
t=0.93440000000000889, A=84.096000000000799, B=0.45517900560517349, AB=1.9709590626248996e-009, 未满足精度要求，返回值仅作参考
t=0.95060000000000922, A=85.554000000000826, B=2.094255315262547e-002, AB=-0.12636660045044271, 未满足精度要求，返回值仅作参考
t=0.95600000000000918, A=86.04000000000083, B=0.21924214499356737, AB=-0.20744208840884626, 未满足精度要求，返回值仅作参考
t=1.0208000000000099, A=91.872000000000895, B=2.7084375215227774e-002, AB=-0.12658892906628622, 未满足精度要求，返回值仅作参考
t=1.0730000000000106, A=96.57000000000096, B=18.643201298211363, AB=7.0910900783409503e-002, 未满足精度要求，返回值仅作参考
t=1.0748000000000106, A=96.732000000000951, B=-0.18505571601464588, AB=1.3235891339475725e-002, 未满足精度要求，返回值仅作参考
t=1.0892000000000108, A=98.028000000000972, B=0.52378102228240464, AB=-0.75649543733381863, 未满足精度要求，返回值仅作参考
t=1.0946000000000109, A=98.514000000000976, B=-0.36850786603778174, AB=-0.14378365430067502, 未满足精度要求，返回值仅作参考
t=1.1594000000000118, A=104.34600000000106, B=0.68701881757929151, AB=-1.5608272756459918, 未满足精度要求，返回值仅作参考
t=1.1684000000000119, A=105.15600000000107, B=2.2634862081548551e-002, AB=-0.21240085279211854, 未满足精度要求，返回值仅作参考
t=1.1702000000000119, A=105.31800000000106, B=-131.97337017437312, AB=-0.22676658336386832, 未满足精度要求，返回值仅作参考
t=1.2134000000000125, A=109.20600000000113, B=-0.39855143896244105, AB=0.36595550982703884, 未满足精度要求，返回值仅作参考
t=1.2296000000000127, A=110.66400000000114, B=-1.112108171122502, AB=4.4630965589931293e-013, 未满足精度要求，返回值仅作参考
t=1.2836000000000134, A=115.52400000000121, B=-1.2812908508775211, AB=-4.4408920985006262e-016, 未满足精度要求，返回值仅作参考
t=1.3016000000000136, A=117.14400000000123, B=-53.25533861126052, AB=7.2997163869104043e-015, 未满足精度要求，返回值仅作参考
t=1.3682000000000145, A=123.13800000000131, B=0.19443836096053715, AB=-9.7590449757639641e-004, 未满足精度要求，返回值仅作参考
t=1.3700000000000145, A=123.3000000000013, B=-0.89667047469271921, AB=1.1043573394329849, 未满足精度要求，返回值仅作参考
t=1.3718000000000146, A=123.46200000000131, B=6.3513193238593608, AB=-0.12505451527879816, 未满足精度要求，返回值仅作参考
t=1.3790000000000147, A=124.11000000000132, B=6.2231361188493777, AB=-0.23863817598755682, 未满足精度要求，返回值仅作参考
t=1.424000000000015, A=128.16000000000136, B=-0.17877662471858569, AB=9.3165938521178604e-003, 未满足精度要求，返回值仅作参考
t=1.4420000000000155, A=129.78000000000139, B=-0.44582445965855844, AB=3.6137759451548845e-014, 未满足精度要求，返回值仅作参考
t=1.4456000000000155, A=130.10400000000141, B=1.7978894937407237e-002, AB=-0.20336766506098625, 未满足精度要求，返回值仅作参考
t=1.5176000000000165, A=136.58400000000148, B=-19.234921536903737, AB=-1.3045120539345589e-014, 未满足精度要求，返回值仅作参考
t=1.5788000000000173, A=142.09200000000155, B=-0.15026132137726742, AB=-0.28323578087700996, 未满足精度要求，返回值仅作参考
t=1.5806000000000173, A=142.25400000000155, B=-0.45099567483060737, AB=2.2204460492503131e-016, 未满足精度要求，返回值仅作参考
t=1.7174000000000194, A=154.56600000000174, B=0.18715643074958366, AB=-4.1211059817981566e-003, 未满足精度要求，返回值仅作参考
t=1.7210000000000192, A=154.89000000000172, B=-0.40729556696552727, AB=-9.5128498133236539e-002, 未满足精度要求，返回值仅作参考
t=1.7696000000000196, A=159.26400000000177, B=-14.047063046372154, AB=-9.6205639518237511, 未满足精度要求，返回值仅作参考
t=1.8578000000000212, A=167.2020000000019, B=0.1606356099221761, AB=-2.646791400684162e-002, 未满足精度要求，返回值仅作参考
t=1.8650000000000211, A=167.8500000000019, B=-1.1617374387243893, AB=-4.8849813083506888e-015, 未满足精度要求，返回值仅作参考
t=1.9280000000000219, A=173.52000000000197, B=0.12750452285190605, AB=-4.5625898014086674e-002, 未满足精度要求，返回值仅作参考
t=1.9298000000000219, A=173.68200000000198, B=12.736580657400062, AB=-8.355504525480778e-002, 未满足精度要求，返回值仅作参考
t=1.9406000000000221, A=174.65400000000199, B=-0.31857494446616447, AB=-9.8693406068930956e-004, 未满足精度要求，返回值仅作参考
t=1.9820000000000226, A=178.38000000000204, B=-0.24243568854155387, AB=4.8079115981049025e-002, 未满足精度要求，返回值仅作参考
t=2., A=180., B=0.37920542909507221, AB=-0.40582069612118654, 未满足精度要求，返回值仅作参考

而x=1时，仅一个输出：

t=1.2386000000000128, A=111.47400000000115, B=13.937513085017597, AB=-1.9365202293231252, 未满足精度要求，返回值仅作参考

rocwoods · 发表于 2010-4-21 22:43:33

17# 闲云孤鹤 没有那么大，4.05G是包含Linux和Windows两个版本的量。windows下完全安装也就4G左右.再说你也可以去掉一些你不要的模块工具箱

zuobinjstu · 发表于 2010-4-26 15:46:38

8# rocwoods rocwoods你好，上面的问题都是我叫同学帮我问的，因为没有这个论坛的账号，刚有了邀请码，注册了这个论坛。我现在又两个难点，还请赐教。1,）我根据您给我的代码，我把它和我的课题结合起来进行了运行，但我知识水平有所欠缺，没有运行成功，在这个代码中是为了验证B和A是等角速传动的，即B的值在90度上下微小波动。代码如下：syms A B ;
r=37;L=140.5;beta=pi/18;
rho=r/2*(1/cos(beta)-1);
theta=asin(rho/L);
f=(r*sin(A)-rho*sin(3*B))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)-...
sin(beta)*sin(3*B)*tan(B)*sin(theta))-(r*cos(A)-rho*cos(beta)*cos(3*B)+...
rho*sin(beta)*sin(theta))*tan(B)*cos(theta);
%将A,B分别用90*t和B(t)替换，为的是好利用符号diff函数来求对B关于t的隐函数F求导
F = subs(f,{'A','B'},{'90*t','B(t)'});
dFt = diff(F,'t');%一阶导数
%将diff(B(t), t)用dBt替换，为的是下一步方便用solve求解diff(B(t), t)的表达式
dFt = subs(dFt,'diff(B(t), t)','dBt');
dBt = solve(dFt,'dBt');%得到B关于t的一阶导数的表达式
B = @(t) fzero(@(B) (r*sin(A)-rho*sin(3*B))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)-...
sin(beta)*sin(3*B)*tan(B)*sin(theta))-(r*cos(A)-rho*cos(beta)*cos(3*B)+...
rho*sin(beta)*sin(theta))*tan(B)*cos(theta) ,90);
%B关于t的函数
eval(['dBt = @(t) ',char(dBt),';' ])
t = 0:0.2:12;
plot(t,arrayfun(@(T) dBt(T),t) )%画dBt关于t的图
2）如何改变此程序中beta值，
当beta=0,，5*pi/180,10*pi/180,15*pi/180,时，dBt（T）的四种曲线绘制在同一张图上。

账号		自动登录	找回密码
密码			注册

账号		自动登录	找回密码
密码			立即注册

隐函数方程的求解[有代码]

相关帖子

评分

评分

本帖子中包含更多资源

评分

评分

本帖子中包含更多资源

评分

本帖子中包含更多资源