Newtons method with two variables

1 次查看（过去 30 天）
Lukas Lehrman
Lukas Lehrman 2018-3-25
编辑： John D'Errico 2018-3-25
So I'm trying to implement Newtons method to find the coordinates of an extremevalue for a function h, whose first and second partial derivates I've named f and g respectively. My code does not return the x and y-coordinates, and I believe the fault lies in my if-statement. Here's the code:
if true
[x,y] = meshgrid (-3:0.1:3,-3:0.1:3);
h = @(x,y) ((x.^3.*y) + ((5*x.^2).*(y.^2)))./(exp(1).^((x.^2) + 3*y.^4));
f = @(x,y) ((3*x.^2.*y)+(10*y.^2.*x)-((10*y.^2).*(x.^3))-(2*x.^4*y))./(exp(1).^((x.^2)+3*y.^4));
%h är funktionen given i uppgiften. f är den partiella derivatan av h med
%avseende på x. Denna omskriving görs på grund av förenklad notation.
g = @(x,y) (((x.^3)+(10*x.^2.*y))-((60*x.^2).*(y.^5))+((12*x.^3).*(y.^4)))./(exp(1).^((x.^2)+3*y.^4));
%h är funktionen given i uppgiften. g är den partiella derivatan av h med
%avseende på y. Denna omskriving görs på grund av förenklad notation.
%contour(x,y,f(x,y)) (contour och mesh skapar graferna i uppgift 2a. Samma
%kod används här som den i uppgift 2a som skapade
%figurerna i labbrapporten.)
%mesh(x,y,f(x,y))
f1 = @(x,y) ((4*x.^5.*y) + ((20*x.^4).*(y.^2))-(14*x.^3.*y)-((50*x.^2).*(y.^2))+(6*x.*y)+(10*y.^2))./(exp(1).^((x.^2)+3*y.^4));
f2 = @(x,y) ((((120*x.^3).*(y.^5))-(120*y.^5.*x))+(((24*x.^3).*(y.^4))-(36*y.^4.*x))+(20*y-(20*x.^2.*y))-2*x.^3+3*x)./(exp(1).^((x.^2)+3*y.^4));
g1 = @(x,y) ((((24*y.^4).*(x.^4))-(2*x.^4))+(((120*y.^5).*(x.^3))-((20*y.^2).*(x.^2))+(3*x.^2)-((36*y.^4).*(x.^2))-(120*y.^5.*x)+(10*y.^2)))./(exp(1).^((x.^2)+3*y.^4));
g2 = @(x,y) (((720*y.^8).*(x.^2))+((144*y.^7).*(x.^3))-(120*y.^5.*x)-((300*y.^4).*(x.^2))-((60*x.^3).*(y.*3))+(20*x.*y))./(exp(1).^((x.^2)+3*y.^4));
x(1) = 1.1;
y(1) = 0.6;
h = 0.0000001;
X = [f(x(i),y(i)) f2(x(i),y(i)); g(x(i),y(i)) g2(x(i),y(i))];
Y = [f1(x(i),y(i)) f(x(i),y(i)); g1(x(i),y(i)) g(x(i),y(i))];
J = [f1(x(i),y(i)) f2(x(i),y(i)); g1(x(i),y(i)) g2(x(i),y(i))];
for i=1
x(i+1) = x(i) - (det(X)/det(J)); %([f(x(i),y(i)) f2(x(i),y(i)); g(x(i),y(i)) g2(x(i),y(i))] / [f1(x(i),y(i)) f2(x(i),y(i));g1(x(i),y(i)) g2(x(i),y(i))]);
x(i) = x(i+1);
y(i+1) = y(i) - (det(Y)/det(J));
y(i) = y(i+1);
if (f(x(i),y(i)) && g(x(i),y(i)) < h )
valX = x(i);
valY = y(i);
return;
end
end
figure()
hold on;
contour(x,y,f(x,y),[0 0], 'b')
contour(x,y,g(x,y),[0 0], 'r')
end
Note: the final four rows simply draw a figure that shows approximate critical points. Comments in swedish. Thanks!

请先登录，再回答此问题。

回答（1 个）

John D'Errico
John D'Errico 2018-3-25
编辑：John D'Errico 2018-3-25
Why in the name of god and little green apples do you need to write your own code?
NEVER write your own code, when code written by professionals is available. This is especially true when you will write poor code.
ezsurf(h,[-3,3])
Rotate the plot around, and it will be obvious the maximum lies near (1,1). It will also be clear there will be multiple solutions. As well, it looks like there might be two solutions with the same maximum, mirror images across the origin.
fminsearch(@(xy) -h(xy(1),xy(2)),[1 1])
ans =
1.062 0.61742
Or, you could have used the optimization toolbox, which will probably give you a few more correct digits for the time required.
format long g
fminunc(@(xy) -h(xy(1),xy(2)),[1 1])
Local minimum found.
Optimization completed because the size of the gradient is less than
the default value of the optimality tolerance.
<stopping criteria details>
ans =
1.06206524800884 0.617438162285868
Next, is it true that the two solutions are the mirrored across the origin?
[xy1,fval] = fminsearch(@(xy) -h(xy(1),xy(2)),[-1 -1])
xy1 =
-1.06203291314256 -0.61741902622501
fval =
-0.604829999447586
[xy1,fval] = fminsearch(@(xy) -h(xy(1),xy(2)),[1 1])
xy1 =
1.06203291314256 0.61741902622501
fval =
-0.604829999447586
Looks like it.
Or, you could use the symbolic TB.
syms x y
H = ((x.^3.*y) + ((5*x.^2).*(y.^2)))./(exp(1).^((x.^2) + 3*y.^4));
Hgrad = gradient(H);
xysol = solve(Hgrad,x,y)
xysol =
struct with fields:
x: [13×1 sym]
y: [13×1 sym]
So 13 solutions, many of which are complex. Which one givex the max?
vpa([xysol.x,xysol.y,h(xysol.x,xysol.y)],8)
ans =
[ 0.90244248, -0.66672465, 0.32319386]
[ -1.0620653, -0.61743817, 0.60483]
[ 1.4133748, -0.14100149, -0.027034788]
[ -1.4133748, 0.14100149, -0.027034788]
[ 1.0620653, 0.61743817, 0.60483]
[ -0.90244248, 0.66672465, 0.32319386]
[ 1.0188998 + 0.076070383i, - 0.025239699 - 0.63521774i, - 0.44731804 - 0.14437493i]
[ 1.0188998 - 0.076070383i, - 0.025239699 + 0.63521774i, - 0.44731804 + 0.14437493i]
[ 10.20644i, -2.0509926i, 1.5918468e23]
[ -10.20644i, 2.0509926i, 1.5918468e23]
[ - 1.0188998 + 0.076070383i, 0.025239699 - 0.63521774i, - 0.44731804 + 0.14437493i]
[ - 1.0188998 - 0.076070383i, 0.025239699 + 0.63521774i, - 0.44731804 - 0.14437493i]
[ 0, 0, 0]
Of the real solutions found, it looks like the two we found before are it.

类别

在 Help Center 和 File Exchange 中查找有关 Calculus 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by