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

huright

function y=main()
    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));
    x0=[60 0];%R,fai
    x=newtons('fRfai',x0,gama,Re,fai0,Au,Rg,p,0.000001)
function [x,fx,xx] = newtons(y,x0,gama,Re,fai0,Au,Rg,p,TolX,MaxIter,varargin)
h = 1e-6; TolFun = eps; EPS = 1e-8;
fx = feval('fRfai',x0,gama,Re,fai0,Au,Rg,p,varargin{:});
Nf = length(fx);
Nx = length(x0);
if Nf ~= Nx, error('Incompatible dimensions of f and x0!'); end
if nargin < 4
    MaxIter = 1000;
end
if nargin < 3
    TolX = EPS;
end
xx(1,:) = x0(:).'
MaxIter = 10;
for k = 1: MaxIter
dx = -jacob('fRfai',xx(k,:),h,gama,Re,fai0,Au,Rg,p,varargin{:})\fx(:);
xx(k + 1,:) = xx(k,:) + dx.';
fx = feval('fRfai',xx(k + 1,:),gama,Re,fai0,Au,Rg,p,varargin{:});
fxn = norm(fx);
if fxn < TolFun | norm(dx) < TolX
    break;
end
end
x = xx(k + 1,:);
if k == MaxIter, fprintf('The best in %d iterations\n',MaxIter), end
function y=fRfai(x,Au,gama,Re,fai0,Rg,p)
    y(1) = (Au - x(1) * cos(x(2))) ^ 2 + (x(1) * sin(x(2)) * cos(gama) - (Re - x(1)) * tan(fai0) * sin(gama)) ^ 2 - Rg ^ 2
    y(2) = (Au + p * tan(gama)) * tan(fai0) * sin(x(2)) - ((Re - x(1)) * tan(fai0) * tan(fai0) -x(1)) * tan(gama) * cos(x(2)) + p - Au * tan(gama)
function g = jacob(f,x,h,gama,Re,fai0,Au,Rg,p,varargin) %Jacobian of f(x)
if nargin < 3
    h = 1e-8;
end
h2 = 2*h;
N = length(x);
x = x(:).';
I = eye(N);
for n = 1:N
g(:,n) = (feval('fRfai',x + I(n,:)*h,gama,Re,fai0,Au,Rg,p,varargin{:})-feval('fRfai',x - I(n,:)*h,gama,Re,fai0,Au,Rg,p,varargin{:}))
end高手们，我用牛顿法编写了程序，怎么解不出来啊？大家帮着看看。

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

xjzuo

稍微修改了一下,看看是不是你想要的结果.
以下程序存为mainprog.m文件,然后在主命令窗口运行之.
%%%==============================%%%
function mainprog
    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));
    x0=[60 0];  %R,fai
    x=mynewtons('fRfai',x0,gama,Re,fai0,Au,Rg,p,0.000001)

%%--------------------------------------------------------------------------------------------------------------%%
function [x,fx,xx] = mynewtons(y,x0,gama,Re,fai0,Au,Rg,p,TolX,MaxIter,varargin)
h = 1e-6; TolFun = eps; EPS = 1e-8;
fx = feval('fRfai',x0,gama,Re,fai0,Au,Rg,p,varargin{:});
Nf = length(fx);
Nx = length(x0);
if Nf ~= Nx, error('Incompatible dimensions of f and x0!'); end
if nargin < 4
    MaxIter = 1000;
end
if nargin < 3
    TolX = EPS;
end
xx(1,:) = x0(:).';
MaxIter = 10;
for k = 1: MaxIter
dx = -myjacob('fRfai',xx(k,:),h,gama,Re,fai0,Au,Rg,p,varargin{:})\fx(:);
xx(k + 1,:) = xx(k,:) + dx.';
fx = feval('fRfai',xx(k + 1,:),gama,Re,fai0,Au,Rg,p,varargin{:});
fxn = norm(fx);
if fxn < TolFun | norm(dx) < TolX
    break;
end
end
x = xx(k + 1,:);
if k == MaxIter, fprintf('The best in %d iterations\n',MaxIter), end

%%-----------------------------------------------------------------------------------------%%
function y=fRfai(x,Au,gama,Re,fai0,Rg,p)
    y(1) = (Au - x(1) * cos(x(2))) ^ 2 + (x(1) * sin(x(2)) * cos(gama) - (Re - x(1)) * tan(fai0) * sin(gama)) ^ 2 - Rg ^ 2;
    y(2) = (Au + p * tan(gama)) * tan(fai0) * sin(x(2)) - ((Re - x(1)) * tan(fai0) * tan(fai0) -x(1)) * tan(gama) * cos(x(2)) + p - Au * tan(gama);

%%------------------------------------------------------------------------------------------------------%%
function g = myjacob(f,x,h,gama,Re,fai0,Au,Rg,p,varargin)  %Jacobian of f(x)
if nargin < 3
    h = 1e-8;
end
h2 = 2*h;
N = length(x);
x = x(:).';
I = eye(N);
for n = 1:N
g(:,n) = (feval('fRfai',x + I(n,:)*h,gama,Re,fai0,Au,Rg,p,varargin{:})-feval('fRfai',x - I(n,:)*h,gama,Re,fai0,Au,Rg,p,varargin{:}))'/h2;
end
%%%==========================================================================%%%

huright

楼上，谢谢。问题解决了。

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