How to make a plot of function in MATLAB

133 Views Asked by At

I have a function ${\dot{\varphi } = \gamma - F(\varphi )}$ (where $F(\varphi)$ - is 2${\pi }$-periodic function) and graph of function $F(\varphi)$.

So I need to create a program that outputs 6 plots of φ(t) for different meaning of $\gamma (\gamma= 0.1, 0.5, 0.95, 1.05, 2, 5)$, and t ∈ [0, 100]

I've got a following F(φ) function: $$F(\phi)=\begin{cases} & {-\frac{\phi}{a}-\frac{\pi}{a}}, &\text{if } {\phi \in [-\pi, -\pi + a]},\\ & -1, &\text{if } {\phi \in [-\pi + a, - a]}, \\ & -\frac{\phi}{a}, &\text{if } {\phi \in [- a, a]}, \\ & 1, &\text{if } {\phi \in [a, \pi- a]},\\ & {-\frac{\phi}{a}+\frac{\pi}{a}},&\text{if } {\phi \in [\pi-a, \pi]}. \end{cases}$$

                 ^
                 |
                 |1   ______
                 |   /|     \
                 |  / |      \
                 | /  |       \
__-π_______-a____|/___|________\π____>
   \        |   /|0    a
    \       |  / |
     \      | /  |
      \     |/   |
       ¯¯¯¯¯¯    |-1

I've tried to do it in MATLAB but I faced with some problems. I don't know what to put in ode45 function (on which section should each part of the graph be displayed and what is the initial value to take). Because the evolution of $\phi(t)$ must be continuous. It's a code for $\gamma= 0.1$

hold on;
df1dt=@(t,f1) 0.1 - f1 - 3.14;
df2dt= @(t,f2)- 1;
df3dt=@(t,f3) 0.1 + f3;
df4dt= @(t,f4)+1;
df5dt=@(t,f5) 0.1- f5 + 3.14;
[T1,Y1] = ode45(df1dt, ...);
[T2,Y2] = ode45(df2dt, ...);
[T3,Y3] = ode45(df3dt, ...);
[T4,Y4] = ode45(df4dt, ...);
[T5,Y5] = ode45(df5dt, ...);
plot(T1,Y1);
plot(T2,Y2);
plot(T3,Y3);
plot(T4,Y4);
plot(T5,Y5);
hold off; 
title('\gamma  = 0.1')
2

There are 2 best solutions below

2
On BEST ANSWER

As a first step, you have to implement the periodicity of $F$ by reducing the argument to a basic period. Your function has odd symmetry, so use $F(-x)=-F(x)$ to reduce the number of cases.

function y = F(x,a)
    n = round(x/(2*pi));
    x = x-n*2*pi; % now in [-pi,pi]
    s = 1; 
    if x<0  x=-x; s=-1; end%if % now in [0,pi]
    y = s*min([1,x/a,(pi-x)/a]);
end%function

Then you can solve your ODE in just one call,

[t,phi] = ode45(@(t,x)gamma-F(x,a), [t0,tf], phi0, options);
plot(t,phi)

The full script that contains the main program and the functions as subroutines in one file can look like this:

function mathse3954790_ballbeam_heun()

  g = 10;
  tau = 0.2;
  J = 1.2;
  m = 5;
  x1_ini = 0.5;
  t0 = 0;
  tf = 10;
  h = 0.05;
  
  x0 = [x1_ini, 0, 0, 0 ];
  t = t0:h:tf;
  x = odeheun(@(t,x)ballbeam(t,x),t,x0); % possibly @ballbeam
  plot(t,x(:,1));

  function dotx = ballbeam(t,x)
      dx2 = x(1)*x(4)^2 - g*sin(x(3));
      dx4 = (-2*m*x(1)*x(2)*x(4)-m*g*x(1)*cos(x(3))+tau)/(m*x(1)^2+J);
      dotx = [ x(2), dx2, x(4), dx4 ];
  end

  function x = odeheun(f,t,x0)
      x = [x0];
      for i=1:(length(t)-1)
          h = t(i+1)-t(i);
          k1 = h*f(t(i),x(i,:));
          k2 = h*f(t(i+1),x(i,:)+k1);
          x(i+1,:) = x(i,:)+0.5*(k1+k2);
      end%for
  end%function

end

There are some tricks that allow to leave out the enclosing function, but it is not very intuitive as the default behavior on loading the script is to call the first function in the script without arguments.

2
On

I have never answered a MATLAB question here. Please ask such type of questions somewhere else.

clear
stp=0.001;
fi=-pi:stp:pi;
a=1;
for s=1:length(fi)
    x=fi(s);
    if x<-pi+a
        f(s)=-(x/a+pi/a);
    elseif x<-a
        f(s)=-1;
    elseif x<a
        f(s)=x/a;
    elseif x<pi-a
        f(s)=1;
    else
        f(s)=-x/a+pi/a;
    end
end
plot(fi,f,'linewidth',1.5)
xlabel('$\phi$','interpreter','latex')
ylabel('$F(\phi)$','interpreter','latex')
title('$a=1$','interpreter','latex')
set(gca,'fontsize',14)
axis([-pi-.5 pi+.5 -1.5 1.5])
grid on

PS: your graph shows a positive slope around zero, but you function doesn't.