Gradient vector field not perpendicular to contour map

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

What's the input data? and what is your purpose?
The model I'm working with is the Kermack-McKendrick SIR model. With some manipulation you get a neat little differential equation that pops out, the orbits of which can be represented by the equation C = I + S - (r/(p*c))*log(S), where r, p, and c are experimental parameters and C is an arbitrary constant. These are the contours that are generated. The gradient uses the same function (f(I,S) = I + S - (r/(p*c))*log(S)), but the subsequent gradient vector field doesn't seem to agree with the contours.
Hi Cole,
could you supply a set of six input parameters that can be used to illustrate the problem?
Give the input directly.....how did you run code and what you expect?
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

0 votes

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