超越方程组的解法（fsolve/牛顿法）

huright

方程组如下：求解fai和R
    N1^2+N2^2=Rg^2
    M1*sin(fai)-M2*cos(fai)+M3=0    %其中：
    m=2;
    r0=26;
    z1=2;
    fai0=pi/6;
    Rf=r0-1.2*m;
    re=60;
    Au=Rf+re;
    p=m*z1/2;
    Ra=r0+m;
    Rg=Rf+0.2*m;
    gama=atan(z1*m/(2*r0));    M1=(Au+p*tan(gama))*tan(fai0);
    M2=((re-R)*tan(fai0)*tan(fai0)-R)*tan(gama);
    M3=p-Au*tan(gama);
    N1=Au-R*cos(fai)
    N2=R*sin(fai)*cos(gama)-(re-R)*tan(fai0)*sin(gama)
先谢谢了

[ 本帖最后由 huright 于 2007-4-22 18:55 编辑 ]

dingd

不限制范围的话有很多组解：
r: 59.6010491005864
fai: -0.00297128614279245


r: -59.7698432177563
fai: -3.07493425931872

r: 106.511057649367
fai: -0.101643383694055

xjzuo

这种问题用fsolve可以解决．
提示：先建立function，再用fsolve调用之．
当x0=[0.5 ; 1]时，结果为：
[phi; R] =
    0.21197005707182
   15.76912924204933

huright

请问怎么解出来的啊？能贴一下程序吗？谢谢

dingd

用1stOpt算的

huright

惭愧啊，俺不知道1stOpt是什么东东，能解释一下吗？

dingd

是个国产独立的优化软件，与Matlab无关，全局优化能力非常强，也非常容易使用。网上搜一下，很多的。

xjzuo

我已经给了提示, 希望huright能动手写写,有问题再来讨论.
这种简单问题看看fsolve的帮助就能解决的.

huright

俺已经用VB编写了程序，但是感觉太复杂，所以想用matlab编；但是用matlab解方程，俺以前没用过，所以求助各位。

笑石头

将你的两个方程组化简
    clear all
    clc
    syms R fai
    m=2;
    r0=26;
    z1=2;
    fai0=pi/6;
    Rf=r0-1.2*m;
    re=60;
    Au=Rf+re;
    p=m*z1/2;
    Ra=r0+m;
    Rg=Rf+0.2*m;
    gama=atan(z1*m/(2*r0));   
    M1=(Au+p*tan(gama))*tan(fai0)   
    M2=((re-R)*tan(fai0)*tan(fai0)-R)*tan(gama)
    M3=p-Au*tan(gama)
    N1=Au-R*cos(fai)
    N2=R*sin(fai)*cos(gama)-(re-R)*tan(fai0)*sin(gama)
解得
M1 = 48.3553


M2 =20/13-4/39*R


M3 = -4.4308


N1 =418/5-R*cos(fai)


N2 =13/170*R*sin(fai)*170^(1/2)-1/510*(60-R)*3^(1/2)*170^(1/2)
将解得M！M1M2M3N1N2带入索求方程组为：
% N1^2+N2^2=Rg^2
(418/5-R*cos(fai))^2+(13/170*R*sin(fai)*170^(1/2)-1/510*(60-R)*3^(1/2)*170^(1/2))^2= 576
%M1*sin(fai)-M2*cos(fai)+M3=0
6805404261977527/140737488355328*sin(fai)-(20/13-4/39*R)*cos(fai)-288/65=0

编写自己的函数：
function F=myfun(x)
% x(1)=R,x(2)=fai;
F=[(418/5-R*cos(fai))^2+(13/170*R*sin(fai)*170^(1/2)-1/510*(60-R)*3^(1/2)*170^(1/2))^2-576;...
     6805404261977527/140737488355328*sin(fai)-(20/13-4/39*R)*cos(fai)-288/65]

用fsolve求解
x0=[0.5 ; 1];
options=optimset( 'MaxIter',30000, 'MaxFunEvals',10000)
[x,fval]=fsolve(@myfun,x0,options)

结果：
x =

   59.6010
   -0.0030


fval =

  1.0e-008 *

    0.8355
   -0.0000

huright

能不能写成一个.m文件啊？还有前面的带入是不是手工带入的阿/〉？fsolve只能调用函数文件名吗？能不能调用方程？？

[ 本帖最后由 huright 于 2007-4-7 15:33 编辑 ]

笑石头

带入不是手工带入
    clear all
    clc
    syms R fai
    m=2;
    r0=26;
    z1=2;
    fai0=pi/6;
    Rf=r0-1.2*m;
    re=60;
    Au=Rf+re;
    p=m*z1/2;
    Ra=r0+m;
    Rg=Rf+0.2*m;
    gama=atan(z1*m/(2*r0));   
    M1=(Au+p*tan(gama))*tan(fai0)   
    M2=((re-R)*tan(fai0)*tan(fai0)-R)*tan(gama)
    M3=p-Au*tan(gama)
    N1=Au-R*cos(fai)
    N2=R*sin(fai)*cos(gama)-(re-R)*tan(fai0)*sin(gama)
f=[N1^2+N2^2-Rg^2;...
    M1*sin(fai)-M2*cos(fai)+M3]
得到：
f =

[ (418/5-R*cos(fai))^2+(13/170*R*sin(fai)*170^(1/2)-1/510*(60-R)*3^(1/2)*170^(1/2))^2-576]
[                6805404261977527/140737488355328*sin(fai)-(20/13-4/39*R)*cos(fai)-288/65]

arcok

好贴，学习学习！

gxc8208

用牛顿-拉弗森法求解

huright

楼上，牛顿-拉弗森法是经常说的牛顿法或者拟牛顿砝码？我查了很多书，也没有找到牛顿-拉弗森法。能给解释一下吗？