0

我即将解出这个方程(-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 
4

1 回答 1

1

这是一个有趣的问题!说明我的无知,我从来没有见过这么扁平的函数!这是“难以实现”的原因是您正在求解的函数非常平坦,这意味着您的自适应步长变得很大,例如~1e23。另外,EPS的公差有点太难了。0.01 相当草率,对于像这个问题这样僵硬的东西,你需要更严格的东西。我投入 1e-100 只是为了好玩,它会像 1 次迭代一样继续迭代,一直迭代到 80!它仍然是一个糟糕的解决方案,但我认为您可以放心您的程序是正确的。我注意到当它过冲太远时,例如 -4,-4,由于 x 和 y 的镜像值,误差正好变为 0。

我的绘图输出

由于问题看起来很棘手,您可能需要考虑一种更稳健的技术来计算自适应步长。!

于 2013-06-05T16:54:12.047 回答