如何求得非线性三角函数的解析解？

bleakhand

有三个方程
第一个，w^2*U0*E0-Bz^2=By0^2;
第二个，w^2*Ud*Ed-Bz^2=Byd^2;
第三个，By0/U0*cot(By0*(b-h))=-Byd/Ud*cot(Byd*h);

要解By0,Byd,Bz三个未知数，结果只需要前三项，如何计算？fsolve好像求不了解析解？

万分感谢！！

ChaChing

个人水平有限, 用了solve试解
Warning: Explicit solution could not be found.
建议试试数值解吧!

bleakhand

谢谢回复，我赋值的话
U0=4*pi*10^(-7);Ud=U0;E0=8.854*10^(-12);Ed=2.56*E0;b=0.01016;h=b/3;
其中w是从1.5e+11到3e+11的矩阵，假设w=1.5e+11:1e+10:3e+11

这样该如何求？我试了下，总说元素不等，即便我用了"."

[ 本帖最后由 ChaChing 于 2010-4-19 11:39 编辑 ]

ChaChing

建议"求助完整格式：出错代码和出错提示":@)

bleakhand

function F=find_beta(B)
w=1.5e+11:1e+10:3e+11;U0=4*pi*10^(-7);Ud=U0;E0=8.854*10^(-12);Ed=2.56*E0;b=0.01016;h=b/3;
F(1)=w.^2.*U0.*E0-B(1).^2-B(2).^2;
F(2)=w.^2.*Ud.*Ed-B(1).^2-B(3).^2;
F(3)=B(2)./U0.*cot(B(2)*(b-h))+B(3)./Ud.*cot(B(3)*h);

fsolve('find_beta',[1 1 1])
???  In an assignment  A(I) = B, the number of elements in B and
I must be the same.

Error in ==> find_beta at 3
F(1)=w.^2.*U0.*E0-B(1).^2-B(2).^2;

Error in ==> fsolve at 254
            fuser = feval(funfcn{3},x,varargin{:});

Caused by:
    Failure in initial user-supplied objective function evaluation. FSOLVE cannot continue.

而且这个初值我确实不知道是什么，[1 1 1]是随便填的。因为解应该也是矩阵，怎么赋值？

ChaChing

抱歉, 这个之前没自己试过, 是来此才学习的! 现在也无法专心试看看!
请先搜索下, 我记得版面有类似的, 如
http://forum.vibunion.com/forum/thread-18886-1-1.html

bleakhand

谢谢。我之前是搜过的，关键词不对，所以没找到这个。谢谢版主。我研究了下那个贴子，可是我的初值还是给不对，总是提示" Optimizer appears to be converging to a minimum that is not a root:Sum of squares of the function values exceeds the square root of options.TolFun. Try again with a new starting point."

初值到底是什么？是不是我猜的这个解大概是多少？

我还有一个简单的问题。

function F=find_omega_c(w)
U0=4*pi*10^(-7);Ud=U0;E0=8.854*10^(-12);Ed=2.56*E0;b=0.01016;h=b/3;
F=sqrt(E0/U0)*cot(w*sqrt(U0*E0)*(b-h))+sqrt(Ed/Ud)*cot(w*sqrt(Ud*Ed)*h);

fsolve('find_omega_c',[1.5e+11])
Optimization terminated: first-order optimality is less than options.TolFun.

ans =

   1.5000e+11

这个w就大概是在e11这个数量级左右，可是为什么不对呢？？这个式子该怎么解？其实，还是初值问题，希望能解答下对于初值的疑惑。谢谢。

[ 本帖最后由 bleakhand 于 2010-2-27 06:55 编辑 ]

bleakhand

请问版主到底有没有办法可以求解三角函数中的所有解？fsolve解不了？

我搜到了这个 http://forum.vibunion.com/forum/viewthread.php?tid=20177
不明白里面的num2str是什么意思，该怎么用？我用matlab没多久，麻烦版主了。

我现在就是要解这个方程
U0=4*pi*10^(-7);Ud=U0;E0=8.854*10^(-12);Ed=2.56*E0;b=0.01016;h=b/3;
F=sqrt(E0/U0)*cot(w*sqrt(U0*E0)*(b-h))+sqrt(Ed/Ud)*cot(w*sqrt(Ud*Ed)*h);

[ 本帖最后由 bleakhand 于 2010-2-28 13:00 编辑 ]

ChaChing

clc; clear
U0=4*pi*10^(-7);Ud=U0;E0=8.854*10^(-12);Ed=2.56*E0;b=0.01016;h=b/3;
syms w
F=sqrt(E0/U0)*cot(w*sqrt(U0*E0)*(b-h))+sqrt(Ed/Ud)*cot(w*sqrt(Ud*Ed)*h);
ezplot(F,[0,2e+11]); grid on
fzero('find_omega_c',[1.5e+11])

ans =
  1.5091e+011

bleakhand

非常感谢！！！！！！！！

我还有个问题，不知道这样太麻烦版主好不好 :P

%求beta的方程组
function F=find_beta(B,w)
U0=4*pi*10^(-7);Ud=U0;E0=8.854*10^(-12);Ed=2.56*E0;b=0.01016;h=b/3;
F=[B(1)^2+B(3)^2-w^2*U0*E0;B(2)^2+B(3)^2-w^2*Ud*Ed;(B(1)/U0)*cot(B(1)*(b-h))+(B(2)/Ud)*cot(B(2)*h)];

%这个是根据版主给我的那个贴子，来画图求解http://forum.vibunion.com/thread-18886-1-1.html

x0=[500;800;0.01];options = optimset('Display','off');w=1.5e+11:1e+10:3e+11;
for i=1:1e+9:length(w)
ww=w(i);B = fsolve(@(B) find_beta(B,ww),x0,options);
B1(i)=B(1);B2(i)=B(2);B3(i)=B(3);
end
plot(w/(2*pi*2.4e+10),B1,'-b',w/(2*pi*2.4e+10),B2,'-r',w/(2*pi*2.4e+10),B3,'-g')

但是最后得到的图不是我想要的

unexpected plot

应该是这样的

wanted

请问哪里出错了？十分感谢ChaChing版主！！！！！！！！！

[ 本帖最后由 bleakhand 于 2010-3-1 13:08 编辑 ]

bleakhand

matlab central有人给我的回复，说fsolve这样解不了

I didn't investigate the dots, but I investigated the function to see if there
was a different approach to solving it. The answer I found was that the
equations are indeed highly non-linear, with multiple poles, and that fsolve()
is not always going to pick up the pole on the same side based upon your
initial guess -- at around the half-way mark, the pole that is picked up for
your initial guess of 500 has moved far enough to the left that there is a
sudden jump and you start picking up the next pole to the right, over in the
800-900 range. Remember, fsolve only finds _an_ answer: it doesn't promise
that the answer is consistent with the ones nearby.

The third equation is the killer, and a change of adding a single extra
decimal to the proposed solution for B(1) can change the calculated value by a
factor of over 100 -- and that's working with 50 decimal places of precision.
With ordinary double precision (16 digits or so), your round off errors in
calculating the third value are going to overwhelm getting correct solutions.
Your B3 plot is likely to be sheer nonsense.


请问版主有什么其他办法可以画吗

ChaChing

LZ的B1/B2/B3根本只有一个值!?
直觉i=1:1e+9:length(w)有问题!

bleakhand

不是啊，500,800,0这个不是初值吗，但是看central别人的回复，这样用fsolve猜初值解不出来

三个方程
B(1)^2+B(3)^2-w^2*U0*E0;
B(2)^2+B(3)^2-w^2*Ud*Ed;
(B(1)/U0)*cot(B(1)*(b-h))+(B(2)/Ud)*cot(B(2)*h);

其实可以变成只含B(3)和变量w的方程。

[ 本帖最后由 bleakhand 于 2010-3-2 11:47 编辑 ]

ChaChing

试试i=1:1:length(w), 不然loop仅有执行一次!

bleakhand

原帖由 ChaChing 于 2010-3-2 11:53 发表
试试i=1:1:length(w), 不然loop仅有执行一次!

好像可以了！！谢谢！！！