[討論] 請大大幫幫忙求近似解
是要把以下程式修改成2 variable並解出
4.*x.^2-20.*x+1/4.*y.^2+8=0
1/2.*x.*y+2.*x-5.*y+8=0
求近似值
但是我要哪裡修改成二元二次的聯立方程？
以下為m檔 小的是非常新的新手還請大大有點耐心教導 謝謝
function steepest(ff)
n=input('輸入方程式數目n=');
ii=0;
for i=1:n
fprintf('輸入初始值x(%d):',i);
x0(i)=input('');
end
tol=input('輸入允許誤差tol=');
maxit=input('輸入最大計算次數maxit=');
fprintf(' k x(1) x(2) x(3) SO\n');
k=0;
% fprintf('%d %f %f %f %f\n',k,x0(1),x0(2),x0(3),s0);
while (1)
y=feval(ff,x0);
s0=y(1).^2+y(2).^2+y(3).^2;%fprintf('s0=%f\n',s0);
fprintf('%d %f %f %f %f\n',k,x0(1),x0(2),x0(3),s0);
t(1)=0;
h(1)=s0;
xr=x0;
xr(1)=xr(1)+0.0001;
yx=feval(ff,xr);
sx=yx(1).^2+yx(2).^2+yx(3).^2;
sxx=(sx-s0)/0.0001;%fprintf('sxx=%f\n',sxx);
xr=x0;
xr(2)=xr(2)+0.0001;
yy=feval(ff,xr);
sy=yy(1).^2+yy(2).^2+yy(3).^2;
syy=(sy-s0)/0.0001;%fprintf('syy=%f\n',syy);
xr=x0;
xr(3)=xr(3)+0.0001;
yz=feval(ff,xr);
sz=yz(1).^2+yz(2).^2+yz(3).^2;
szz=(sz-s0)/0.0001;%fprintf('szz=%f\n',szz);
M=sqrt(sxx.^2+syy.^2+szz.^2);%fprintf('M=%f\n',M);
if abs(M)<=0.000001
%fprintf('ans=%f
break;
end
ir=1;iim=0;
t(3)=1.0;
%while(ir>=1)
xt(1)=x0(1)-t(3).*sxx/M;
xt(2)=x0(2)-t(3).*syy/M;
xt(3)=x0(3)-t(3).*szz/M;
y=feval(ff,xt);
h(3)=y(1).^2+y(2).^2+y(3).^2;%fprintf('h(1)=%f,h(3)=%f\n',h(1),h(3));
while(ir>=1)
ir=0;
if abs(h(3))>=abs(h(1))
t(3)=t(3)/2;
if t(3)<0.00001
%fprintf('t(3)is too small...\n');iim=1;
break;
end
ir=1;
end
xt(1)=x0(1)-t(3).*sxx/M;
xt(2)=x0(2)-t(3).*syy/M;
xt(3)=x0(3)-t(3).*szz/M;
y=feval(ff,xt);
h(3)=y(1).^2+y(2).^2+y(3).^2;%fprintf('h(1)=%f,h(3)=%f\n',h(1),h(3));
end
if iim==1
break;
end
t(2)=t(3)/2;
h(1)=s0;
xt(1)=x0(1)-t(2).*sxx/M;
xt(2)=x0(2)-t(2).*syy/M;
xt(3)=x0(3)-t(2).*szz/M;
y=feval(ff,xt);
h(2)=y(1).^2+y(2).^2+y(3).^2;
p1=(h(2)-h(1))/(t(2)-t(1));
p2=(h(3)-h(2))/(t(3)-t(2));
p3=(p2-p1)/(t(3)-t(1));
t(4)=0.5.*(t(1)+t(2)-p1/p3);
xt(1)=x0(1)-t(4).*sxx/M;
xt(2)=x0(2)-t(4).*syy/M;
xt(3)=x0(3)-t(4).*szz/M;
y=feval(ff,xt);
h(4)=y(1).^2+y(2).^2+y(3).^2;%fprintf('h(1)=%f,h(2)=%f,h(3)=%f,h(4)=%f\n',h(1),h(2),h(3),h(4));
min=h(1);%fprintf('t(1)=%f,t(2)=%f,t(3)=%f,t(4)=%f\n',t(1),t(2),t(3),t(4));
kk=1;
for j=1:4
if h(j)<min
min=h(j);
kk=j;
end
end
tt=t(kk);hh=h(kk);%fprintf('tt=%f\n',tt);
xt(1)=x0(1)-tt.*sxx/M;
xt(2)=x0(2)-tt.*syy/M;
xt(3)=x0(3)-tt.*szz/M;
%fprintf('hh=%f,xt(1)=%f,xt(2)=%f,xt(3)=%f,so=%f\n',hh,xt(1),xt(2),xt(3),so);
if abs(hh-h(1))<0.0001
break;
end
k=k+1;
if k>=maxit
break;
end
x0=xt;
%fprintf('%d %f %f %f %f\n',k,x0(1),x0(2),x0(3),s0);
end
--
※ 發信站: 批踢踢實業坊(ptt.cc), 來自: 123.192.77.189
※ 文章網址: http://www.ptt.cc/bbs/MATLAB/M.1402500511.A.1C1.html
MATLAB 近期熱門文章
PTT數位生活區 即時熱門文章
