1

I've got the following code in matlab:

gamma=1.01;
[t,phi] = ode45(@(t,x)(gamma-F(x,pi/6, 0.5)), [0,100], 0);
plot(t,phi)
hold off; 
title('\gamma = 1.01')

function out = F(p, a, b)
phi = mod(p,2*pi);
out = (0      <= phi & phi < a     ).*(phi/a) ...
    + (a      <= phi & phi < pi/2  ).*((phi*(b-1)/((pi/2)-a))-(a*(b-1))/((pi/2)-a)+1)...
    + (pi/2   <= phi & phi < pi-a  ).*((phi*(1-b)/((pi/2)-a))-((pi-a)*(1-b))/((pi/2)-a)+1)...
    + (pi-a   <= phi & phi < pi+a  ).*(-phi/a + pi/a)...
    + (pi+a   <= phi & phi < 3*(pi/2)).*((phi*(1-b)/((pi/2)-a))-(((pi+a)*(1-b))/((pi/2)-a))-1)...
    + (3*pi/2 <= phi & phi < 2*pi-a).*((phi*(b-1)/((pi/2)-a))-(((2*pi-a)*(b-1))/((pi/2)-a))-1)...
    + (2*pi-a <= phi & phi < 2*pi  ).*(phi/a-2*pi/a);
end

And I have a question. Could you please say how to numerically integrate (I suppose it's called that way) this function in MATLAB?

So from this enter image description here Get this enter image description here

Jane
  • 21
  • 3

1 Answers1

0

For numerical integration you should probably use simple trapz: https://www.mathworks.com/help/matlab/ref/trapz.html. You can also use integral https://www.mathworks.com/help/matlab/ref/integral.html.

So, for your code, you need simple

area=trapz(t,phi);

This is the whole integral of the function - total area underneath the curve. If you want to show how integral grows with t, you need to use cumtrapz function which does cumulative integration:

area = cumtrapz(t,phi);
plot(t, area)
% This is essentially equivalent to doing `trapz` in a loop with less hassle (for... area(i)=trapz(t(1:i),phi(1:i));...)

However, your red line does not show integral (integral would be area underneath that red or blue line). It seems you want approximation of the function at few points, and draw these points connected? For that, you can try searching for points with local maxima/minima of second derivative, which should correspond to desired points.

To get derivative, use diff: https://www.mathworks.com/help/matlab/ref/diff.html

So

der1 = diff(phi);
der2 = diff(der1);
isLocalMax = (der2(2:end-1) > der2(1:end-2)) & (der2(2:end-1) > der2(3:end));
isLocalMin % equivalent...
plot(t(isLocalMax|isLocalMin), phi(isLocalMax|isLocalMin))
Zizy Archer
  • 1,392
  • 7
  • 11
  • @ZiziArcher thank you a lot! I wanted to ask about you using `trapz`. I added it into code but the plot is the same. Should I change something else? – Jane Oct 29 '21 at 07:24
  • @Jane Integral as given by `trapz` is the whole area (= a single number), I now added solution with cumulative integration `cumtrapz` if you want to see how area grows over time (`t`). – Zizy Archer Oct 29 '21 at 07:40