Evaluate function over a mesh grid (without for loops)

Is there a way to evaluate a function that takes an [x, y] vector as input over a grid of points?
Here is the function I want to evaluate:
Jfun = @(u) (17*u(1)^2)/2 - 14*u(1)*u(2) - 40*u(1) + 19*u(2)^2 - 20*u(2)
The reason I defined the function this way is so that I can find the minimum thus:
umin = fminsearch(Jfun,[0,0])
Now I want to make the contour plot.
Something like this perhaps:
x = linspace(-1,6);
y = linspace(-3,5);
[X,Y] = meshgrid(x,y);
for i=1:numel(x)
for j=1:numel(y)
Z(i,j)=Jfun([X(j,i) Y(j,i)]);
end
end
contourf(X,Y,Z,10)
But is there a way to compute Z without having to use nested for loops?
I tried this:
Z = Jfun(cat(3,X,Y));
but it doesn't work. It simply returns 173.5000 which is Jfun([-1 -1]).

 採用された回答

Rik
Rik 2020 年 10 月 22 日

0 投票

You have to vectorize the function and split it into two steps:
Jfun = @(x,y) (17*x.^2)/2 - 14*x.*y - 40*x + 19*y.^2 - 20*y;
umin = fminsearch(@(u) Jfun(u(1),u(2)),[0,0])
x = linspace(-1,6);
y = linspace(-3,5);
[X,Y]=ndgrid(x,y);
Z=Jfun(X,Y);

3 件のコメント

Bill Tubbs
Bill Tubbs 2020 年 10 月 22 日
Thanks, this works.
Bill Tubbs
Bill Tubbs 2020 年 10 月 22 日
Unfortunately the execution of the fminsearch is about 40% slower with this approach but the computation of Z is very fast.
Rik
Rik 2020 年 10 月 22 日
You can remove the overhead of the outer layer of anonymous function if that 40% is important for you. That does mean you will have two copies of virtually the same function, but that is a small price to pay.

サインインしてコメントする。

その他の回答 (1 件)

madhan ravi
madhan ravi 2020 年 10 月 22 日

0 投票

v = num2cell([X(:), Y(:)], 2);
Z = reshape(cellfun(Jfun, v), size(X));

1 件のコメント

Bill Tubbs
Bill Tubbs 2020 年 10 月 22 日
Thanks. This works but is not very efficient according to timeit().

サインインしてコメントする。

カテゴリ

ヘルプ センター および File ExchangeLoops and Conditional Statements についてさらに検索

製品

リリース

R2019b

質問済み:

2020 年 10 月 22 日

コメント済み:

Rik
2020 年 10 月 22 日

Community Treasure Hunt

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

Start Hunting!

Translated by