高手们，我用牛顿法编写了程序，怎么解不出来啊？大家帮着看看。
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=;%R,fai
x=newtons('fRfai',x0,gama,Re,fai0,Au,Rg,p,0.000001)
function = 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 编辑 ] 稍微修改了一下,看看是不是你想要的结果.
以下程序存为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=;%R,fai
x=mynewtons('fRfai',x0,gama,Re,fai0,Au,Rg,p,0.000001)
%%--------------------------------------------------------------------------------------------------------------%%
function = 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 于 2007-4-22 21:09 编辑 ]
页:
1
[2]