How to plot a function which is defined on different subintervals

I am trying to plot the function ,

But I don't know how to write the code for the definition of the function f which is given on different subintervals.

[Jan]

At the points i, f has two different values. Shouldn't one of the intervals exclude i?

[Cris19]

@JanThank you for the comment. Indeed, f had two different values at i and I have corrected the function. Now it should be continuous.

[KALYAN ACHARJYA]

MATLAB Onramp

[Cris19]

@KALYAN ACHARJYA Thank you, but didn't help.

Answer 1

[Walter Roberson]

There are several methods available. The most straight forward is to write a function that loops over the inputs, testing each one to decide what the result should be.

function y = f(X)
   y = zeros(size(X));
   for K = 1 : numel(X)
       x = X(K);
       if x < 2
          y(K) = 0;
       elseif x <= 5
           y(K) = x.^2;
       elseif x <= 8
           y(K) = x-x.^3;
       else
           y(K) = 0;
       end
end

With regards to those intervals you need, think about floor(x+1/2)

[Cris19]

@Walter Roberson Thank you, but I still don't know how to write the code for f(x) on [2,500]. After defining and plotting the function f(x) on the interval [2,500], I also want to plot on [2,500] the solution of the second order ODE y''+f(x)y'+y=0 with the initial data y(2)=0.1, y'(2)=0.1. This is why I think we should know the function f(x) on the whole interval [2,500].

[Walter Roberson]

The function is defined in sections. You would have a difficult time writing it as a formula unless you use piecewise or an infinite sum of heaviside.

[Cris19]

How exactly can I use the piecewise syntax? Would you please write the code? I am not too familiarized with Matlab.

[Cris19]

I tried to use this code:

for i= 2:500
Bounds_1 = [i - (1/i^2),i];
Bounds_2 = [i,i + (1/i^2)];
Bounds = [Bounds_1 Bounds_2];
Min = min(Bounds);
Max = max(Bounds);
f1 = @(x) (i^2.*x - i^3 + 1).*(Bounds_1(1) < x & x <= Bounds_1(2));
f2 = @(x) (-i^2.*x + i^3 + 1).*(Bounds_2(1) < x & x <= Bounds_2(2));
f = @(x) f1(x) + f2(x); 
end
fplot(f);
xlim([Min-2 Max+2]);
ylim([0 1.1]);

But the only plotted 'spike' is the last one, for i=500. How can I keep all spikes in the expression of f?

[Walter Roberson]

for i= 2:500

i is a positive integer

Bounds_1 = [i - (1/i^2),i];

With i being real valued, 1/i^2 is real valued and positive (or infinite and positive), and subtracting a positive amount from a real value always gives you a result that is less than the original amount.

Bounds_2 = [i,i + (1/i^2)];

With i being real-valued, 1/i^2 is real valued and positive (or infinite and positive), and adding a positive amount to a real value always gives you a result that is greater than the original amount

Bounds = [Bounds_1 Bounds_2];

That would be a list [i - (1/i^2),i,i,i + (1/i^2)]

Min = min(Bounds);

As demonstrated above, i-1/i^2 is always less than i, and i+1/i^2 is always greater than i, so we can quickly see that the minimum value must be i-1/i^2

Max = max(Bounds);

and the maximum value must be 1+1/i^2 . So there is no point in creating the Bounds list and taking the min and max: we already know which will be the min and which will be the max

end

You do not use Min and Max inside the loop, so you overwrite the value each time. Therefore at the end of the loop, the value will be whatever was calculated on the last iteration, the one for i=500

fplot(f);
xlim([Min-2 Max+2]);

and you create your xlim according to only that iteration.

f1 = @(x) (i^2.*x - i^3 + 1).*(Bounds_1(1) < x & x <= Bounds_1(2));
f2 = @(x) (-i^2.*x + i^3 + 1).*(Bounds_2(1) < x & x <= Bounds_2(2));
f = @(x) f1(x) + f2(x); 
end

You overwrite all of f1, f2, and f each iteration, so in that code there is no point doing any iteration except the last one.

[Walter Roberson]

You could try this:

f = @(x) zeros(size(x));
for i= 2:500
    Bounds_1 = [i - (1/i^2),i];
    Bounds_2 = [i,i + (1/i^2)];
    f1 = @(x) (i^2.*x - i^3 + 1).*(Bounds_1(1) < x & x <= Bounds_1(2));
    f2 = @(x) (-i^2.*x + i^3 + 1).*(Bounds_2(1) < x & x <= Bounds_2(2));
    f = @(x) f(x) + f1(x) + f2(x); 
end
fplot(f, [1 6.5]); ylim([-0.1 1.1])

It is difficult to plot much further than this, as the spikes get so narrow that fplot misses them. By 500 the spikes are only 1/2500 wide, and assuming you want at least 5 points around 500, that implies you would have to use a sampling frequency no less than 125000 Hz, times your range of at least 501 (you need to include 1 and 501 exactly in the sample so that the fourier transform does not accidentally think that the pattern repeats indefinitely when instead it stops at 1 and 501.)

[Cris19]

Thank you so very much for your thorough comments! The code works very well!

[Cris19]

One more question, if possible... I try to solve to second order ODE y''+f(x)y'+y=0, whre f(x) is the function above. I am interested to plot on [0,10] the solution of this ODE. for the initial data y(0)=y'(0)=0.1. I used the function (where f(x) is defined with the code you wrote)

function dwdt=susp(t,w)
f = @(x) zeros(size(x));
for i= 2:10
Bounds_1 = [i - (1/i^2),i];
Bounds_2 = [i,i + (1/i^2)];
f1 = @(x) (i^2.*x - i^3 + 1).*(Bounds_1(1) < x & x <= Bounds_1(2));
f2 = @(x) (-i^2.*x + i^3 + 1).*(Bounds_2(1) < x & x <= Bounds_2(2));
f = @(x) f(x) + f1(x) + f2(x);
end
dwdt=zeros(2,1);
dwdt(1)=w(2);
dwdt(2)=-f*w(2)-w(1);
end 

and the commands

tspan = [0 10];
z0=[0.1 0.1];
[t,z] = ode45(@(t,z) susp(t,z), tspan, z0);
figure
plot(t,z(:,1),'r');
legend('x(t)');

But I get several errors...

Unary operator '-' is not supported for operand of type 'function_handle'.

Error in susp (line 12)
dwdt(2)=-f*w(2)-w(1);

Error in @(t,z)susp(t,z)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

How could I fix them?

[Walter Roberson]

It is not possible to use the ode*() functions with discontinuous functions (or rather, if it works at all it will give the wrong answer.)

The exception to that is if the ode*() call has continuous second derivatives for the entire call then it is allowed. So you would have to break up the time into segments: [2, 2+1/4], [2+1/4, 3-1/9], [3-1/9, 3], [3,3+1/9], [3+1/9, 4-1/16], [4-1/16, 4], [4,4+1/6] and so on.

[Cris19]

Thank you. I wonder if there is an alternative solution. I am interested to find the asymptotic behavior at +infinity of the solution of that ODE. So is this why I think it is interesting to see the plotting of the solution on intervals larger and larger, such as [0,10], [0,500] etc.

I really don't know how to code this...

[Walter Roberson]

At asymptopic behaviour is

syms n f(x)
D = symsum(dirac(x-n), n, 1, 200)
df = diff(f)
d2f = diff(df)
eqn = d2f + D*df + f == 0
dsolve([eqn, f(0)==0])

but in practice you will not get any solution. Even if you reduce it down to dirac at one particular integer you are not going to get a solution.