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

zuobinjstu

20# zuobinjstu A=90，B也应该是90左右波动，我经行了运动仿真，B根据beta角度不同的仿真图形如附件。我想用公式理论验证下

rocwoods

1. 
   
2. B = @(t) fzero(@(B) (r*sin(A)-rho*sin(3*B))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)-...
   
3. sin(beta)*sin(3*B)*tan(B)*sin(theta))-(r*cos(A)-rho*cos(beta)*cos(3*B)+...
   
4. rho*sin(beta)*sin(theta))*tan(B)*cos(theta) ,90);
   

复制代码

改为

1. 
   
2. B = @(t) fzero(@(B) (r*sin(90*t)-rho*sin(3*B))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)-...
   
3. sin(beta)*sin(3*B)*tan(B)*sin(theta))-(r*cos(90*t)-rho*cos(beta)*cos(3*B)+...
   
4. rho*sin(beta)*sin(theta))*tan(B)*cos(theta) ,90);
   

复制代码

就可以了。关于如何改变此程序中beta值，整段代码做成函数形式，接受bete输入，out = arrayfun(@(T) dBt(T),t); out作为对应beta值的输出，这样对应不同的beta得到不同的out，然后画图就OK了。

zuobinjstu

22# rocwoods 谢谢你的多次帮忙，让我受益匪浅。也感谢这个论坛有这么多的热心朋友的帮助。你也一定很忙吧！那么晚了还回复我，我很感谢！rocwoods如果可以，可以加我为好久嘛？我的邮箱：[email=zuobinjstu@163.com, QQ:354637093]zuobinjstu@163.com, QQ:354637093[/email]

zuobinjstu

22# rocwoods
我运行出来图像了，可让我纳闷的是，这个图像和理论所要的图像完全不一样，不知道怎么回事。我的导师曾经将那个公式用手算化简出来然后汇出的图像应该在90左右，（见20楼附件）而我运行出来图像完全不一样，我在想难道是程序哪里有偏差还是matlab运算图像的算法不一样。还请你帮忙看下，谢谢！
代码如下是为了算这个：

当(r*sin(90*t)-rho*sin(3*B)*cos(theta))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)*cos(theta)-...

sin(beta)*sin(3*B)*tan(B)*sin(theta))-(r*cos(90*t)-rho*cos(beta)*cos(3*B)+...

rho*sin(beta)*sin(theta))*tan(B)*cos(theta)=0时， B的一阶导数关于t的图像（应该是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(theta))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)*cos(theta)-...

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

F = subs(f,{'A','B'},{'90*t','B(t)'});

dFt = diff(F,'t');

dFt = subs(dFt,'diff(B(t), t)','dBt');

dBt = solve(dFt,'dBt');

B = @(t) fzero(@(B) (r*sin(90*t)-rho*sin(3*B)*cos(theta))*(cos(beta)-...

sin(beta)*sin(theta)*cos(3*B)*cos(theta)-...

sin(beta)*sin(3*B)*tan(B)*sin(theta))-(r*cos(90*t)-rho*cos(beta)*cos(3*B)+...

rho*sin(beta)*sin(theta))*tan(B)*cos(theta) ,90);

eval(['dBt = @(t) ',char(dBt),';' ])

t = 0:0.5:12;

plot(t,arrayfun(@(T) dBt(T),t) )

rocwoods

你这个问题原因在于方程是高频振荡方程，有多个解。而fzero求解只能求出一个解，fzero求解速度非常快，但是不能保证求得的解在初始点附近，因此出现了有些解远离90的情况。只能退而求其次用fsolve替换fzero，fsolve本来是用来求多元非线性方程的解的，给定初始值，一般求出的解不会脱离初值太远，但是速度比较慢。修改后代码如下：

1. 
   
2. syms A B ;
   
3. r=37;L=140.5;beta=pi/18;
   
4. rho=r/2*(1/cos(beta)-1);
   
5. theta=asin(rho/L);
   
6. f=(r*sin(A)-rho*sin(3*B)*cos(theta))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)*cos(theta)-...
   
7. sin(beta)*sin(3*B)*tan(B)*sin(theta))-(r*cos(A)-rho*cos(beta)*cos(3*B)+...
   
8. rho*sin(beta)*sin(theta))*tan(B)*cos(theta);
   
9. F = subs(f,{'A','B'},{'90*t','B(t)'});
   
10. dFt = diff(F,'t');
    
11. dFt = subs(dFt,'diff(B(t), t)','dBt');
    
12. dBt = solve(dFt,'dBt');
    
13. 
    
14. options = optimset('TolFun',1e-12,'display','off');
    
15. B = @(t) fsolve(@(B) (r*sin(90*t)-rho*sin(3*B)*cos(theta)).*(cos(beta)-...
    
16. sin(beta)*sin(theta)*cos(3*B)*cos(theta)-...
    
17. sin(beta)*sin(3*B).*tan(B)*sin(theta))-(r*cos(90*t)-rho*cos(beta)*cos(3*B)+...
    
18. rho*sin(beta)*sin(theta)).*tan(B)*cos(theta) ,90,options);
    
19. 
    
20. eval(['dBt = @(t) ',char(dBt),';' ])
    
21. t = 0:0.5:12;
    
22. plot(t,arrayfun(@(T) dBt(T),t) )
    

复制代码

如果还是不符合预期，试着从你专业角度入手，看看方程可以怎样简化或者变形，使得方程的解能够限制在可控的范围内。

zuobinjstu

25# rocwoods 恩，这个图像是比fzero好多了，我需要的图像是在在90之间微微波动，很显然还需要在处理。那个f函数是这个，我们输入到matlab里的表达式：

• f=(r*sin(A)-rho*sin(3*B)*cos(theta))*(cos(beta)-sin(beta)*sin(theta)*cos(3*B)*cos(theta)-...
• 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)与下图是一样的吧

zuobinjstu

25# rocwoods 好像很难根据目前的表达式求出以下附件的图形。90左右微微波动