Gradient vector field not perpendicular to contour map

6 vues (au cours des 30 derniers jours)
Cole Butler
Cole Butler le 10 Oct 2017
Commenté : Cole Butler le 16 Oct 2017
I have the following code:
function [ ] = kermack( I0,r,p,c,N,numSteps )
I(1)=I0;
S(1)=N-I0;
for t = 1:numSteps % construct I and S vectors
I(t+1) = I(t) + (p*c*I(t)*S(t) - r*I(t));
S(t+1) = S(t) - p*c*I(t)*S(t);
end
% create contour
x = min(I):0.2:max(I);
y = min(S):0.2:max(S);
[xx,yy] = meshgrid(x,y);
zz = xx + yy - (r/(p*c))*log(yy); % f(I,S)
[DX,DY] = gradient(zz);
figure % new window
contour(xx,yy,zz)
hold on % superimpose graphs
scale=10; % to make the arrows visible
quiver(xx,yy,DX,DY,scale) % vector field of f(I,S)
xlim([0 60]); % set bounds on the window
ylim([0,80]);
xlabel('I(t) (infected/infectious population)');
ylabel('S(t) (susceptible population)');
hold off
Now, as made clear by the title of my question, the generated gradient vector field doesn't seem to be perpendicular with the superimposed contour mapping. How can this be? What am I doing wrong?
  5 commentaires
KSSV
KSSV le 12 Oct 2017
Give the input directly.....how did you run code and what you expect?
Cole Butler
Cole Butler le 16 Oct 2017
In response to the previous comments, I have attached an example of the code. This is the output of kermack() for p=0.05, c=0.2, r=0.2, numSteps = 50, N=100, I0=10. Clearly, the solid blue lines are the contour lines. When you zoom in, you can see that the gradient field is not perpendicular to the contour lines.
<<
>>

Connectez-vous pour commenter.

Réponse acceptée

David Goodmanson
David Goodmanson le 16 Oct 2017
Hi Cole,
I believe this is an artifact of the zoom process. If you zoom in by drawing a rectangle, the axes get scaled differently and right angles no longer appear to be right angles. If you add
axis equal
at the end of the code and then zoom in simply by repeated clicking on a point on a contour line (which does 2x on each axis), the results appear to be good.

Plus de réponses (0)

Community Treasure Hunt

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

Start Hunting!

Translated by