我即将解出这个方程(-cos(x).*cos(y).*exp(-((x-pi).^2+(y-pi).^2)=0
。我的代码似乎适用于其他更简单的方程,但这个很难实现。有什么建议么?
clc
close
clear
[x,y] = meshgrid([-20:0.1:20],[-20:0.1:20]);
fn = (-cos(x).*cos(y).*exp(-((x-pi).^2+(y-pi).^2)));
f=@(x,y) (-cos(x).*cos(y).*exp(-((x-pi).^2+(y-pi).^2)));
h=0.001;
dfx=@(x,y) (f(x-2*h,y)-8*f(x-h,y)+8*f(x+h,y)-f(x+2*h,y))/(12*h);
dfy=@(x,y) (f(x,y-2*h)-8*f(x,y-h)+8*f(x,y+h)-f(x,y+2*h))/(12*h);
d2fx=@(x,y) (dfx(x-2*h,y)-8*dfx(x-h,y)+8*dfx(x+h,y)-dfx(x+2*h,y))/(12*h);
d2fy=@(x,y) (dfy(x,y-2*h)-8*dfy(x,y-h)+8*dfy(x,y+h)-dfy(x,y+2*h))/(12*h);
dfxfy=@(x,y) (dfx(x,y-2*h)-8*dfx(x,y-h)+8*dfx(x,y+h)-dfx(x,y+2*h))/(12*h);
dfyfx=@(x,y) (dfy(x-2*h,y)-8*dfy(x-h,y)+8*dfy(x+h,y)-dfy(x+2*h,y))/(12*h);
subplot(2,1,1)
mesh(x,y,f(x,y));
subplot(2,1,2)
contour(x,y,f(x,y))
hold on
l= 1;
dx(1) = -2;
dy(1) = -2;
k = 1;
licznik = 1;
while( 1==1 )
xg = [dfx(dx(l),dy(l))];
yg= [dfy(dx(l),dy(l))];
a=([xg;yg]'*[xg;yg])...
/([xg;yg]'*[d2fx(dx(l),dy(l)),dfxfy(dx(l),dy(l));dfyfx(dx(l),dy(l)),d2fy(dx(l),dy(l))]*[xg;yg]);
wyn = [dx(l);dy(l)]-a*[xg;yg];
dx(l+1)=wyn(1);
dy(l+1)=wyn(2);
if(abs(f(dx(l+1),dy(l+1)) -f(dx(l),dy(l)))<0.01 || l >=100)
break;
end
l=l+1;
end
plot(dx(1),dy(1),'r*')
text(dx(1),dy(1),'START')
plot(dx(2:length(dx)-1),dy(2:length(dy)-1),'g.')
plot(dx(end),dy(end),'r*')
text(dx(end),dy(end),'STOP')
for i=1:length(dx)-1
line([dx(i),dx(i+1)],[dy(i),dy(i+1)],'color','g')
end
for i=1:length(dx)
contour(x,y,f(x,y),[f(dx(i),dy(i)),f(dx(i),dy(i))])
end