Using ode45 and plot

I am working on 1. I am having trouble coding the function so that I can use the ode45 function. This is what my code looks like :

function xp=F(t,x)
xp=zeros(2,1); % since output must be a column vector
xp(2)=-(x(1)+(0*(x(1)^3)));
xp(3)=-(x(1)+(0.1*(x(1)^3)));
xp(4)=-(x(1)+(0.2*(x(1)^3)));
xp(5)=-(x(1)+(0.4*(x(1)^3)));
xp(6)=-(x(1)+(0.6*(x(1)^3)));
xp(7)=-(x(1)+(0.8*(x(1)^3)));
xp(8)=-(x(1)+(1.0*(x(1)^3)));

then in the command window I am inputting :

[t,x]=ode45('F',[0,20],[0,1]); plot(t,x(:,1))

I know I probably messed up the code . Any ideas?

Answer 1

[Star Strider]

Don’t do everything at once!

Otherwise, you’re almost there! This is how I would do it:

F = @(t,x,e) [x(2); -(x(1)+(e.*(x(1).^3)))];        % Spring ODE

ev = 0:0.2:1.0;                                     % ‘epsilon’ Vector
ts = linspace(0,10, 50);                            % Time Vector

for k1 = 1:length(ev)
    [~, x{k1}]=ode45(@(t,x) F(t,x,ev(k1)), ts, [0,1]);
end

figure(1)
for k1 = 1:length(ev)
    subplot(length(ev), 1, k1)
    plot(ts,x{k1})
    legend(sprintf('\\epsilon = %.1f', ev(k1)))
    grid
    axis([xlim    -2  +2])
end
xlabel('Time')

I stored ‘x’ as a cell array because it’s easier than storing it as a multi-dimensional array. The first variable, ‘x(:,1)’ is blue and ‘x(:,2)’ is red in each plot. I used subplots because it’s easier to compare the plots that way. I also specified the time argument, ‘ts’ as a vector so that all the integrated values would be the same size. These choices are for convenience, and not required.

[Sam]

Thank you for the help! So looking at the graphs would the amplitude just be 1? Also, I graphed epsilon values going to 50 and the graph was basically a straight line on 0. Would this mean that as epsilon --> infinity the amplitude goes to 0?

Thanks!

[Star Strider]

My pleasure!

I got similar results when I (just now) varied epsilon from 0 to 50 in steps of 10. The value of epsilon affects the frequency, from 1/6.25 for epsilon = 0 to about 1/2.25 for epsilon = 50, while the amplitudes of ‘x(:,1)’ decreased from ±1 to ±0.42 as epsilon increased, while ‘x(:,2)’ remained stable at about ±1.

So you are correct — the amplitude of ‘x(:,1)’ varies inversely with epsilon, and the frequency varies proportionally with epsilon.

[Sam]

I am trying to find the value of t when the graph first hits the equilibrium(0) I have been using the data cursor on the plot, but it is not precise enough because I am getting the same values for when epsilon = 0.2 and 0.4 to be 2.857 and I know this shouldnt be the case, it should be decreasing. Is there a function I can use so that it prints these exact values?

thanks!

[Star Strider]

Here is the way I usually use to detect zero-crossings:

t = 10:12;
x = [2 0 -2]';
zx = t(x .* circshift(x,[0 1]) <= 0)

The ‘zx’ assignment here gives the ‘t’ value for the zero-crossing. The ‘x’ vector does not have to include zero (it did here for illustration purposes) but the ‘x’ variable must change signs.

This will give you the approximate location of the zero-crossing. To get a more precise value, use an interpolation routine. For a detailed discussion of the techinque, see Curve Fitting to a Sinusoidal Function.

Note: Column vectors and row vectors have to be matched to the circshift second argument.

———————————————————————————————————

EDIT — This version of my original code will find and plot the initial zero-crossings for x(:,1). (In the ‘ixrng’ assignment, using the third value of ‘zxix’ was necessary for the code to avoid the initial value as the first zero.)

F = @(t,x,e) [x(2); -(x(1)+(e.*(x(1).^3)))];                % Spring ODE

ev = 0:0.2:1.0;                                             % ‘epsilon’ Vector
ev = [0:1:5]*10;
ts = linspace(0,6.5, 50);                                   % Time Vector

for k1 = 1:length(ev)
    [~, x{k1}]=ode45(@(t,x) F(t,x,ev(k1)), ts, [0,1]);
end

figure(1)
for k1 = 1:length(ev)
    subplot(length(ev), 1, k1)
    plot(ts,x{k1})
    hold on
    legend(sprintf('\\epsilon = %.1f', ev(k1)))
    grid
    axis([xlim    -2  +2])
    xv = x{k1}(:,1);
    zxix = find((xv.*circshift(xv, [1 0])) <= 0);           % Approximate Zero-Crossing Indices
    ixrng = max(1,zxix(3)-1):min(zxix(3)+1,length(ts));     % Index Range
    tzx{k1} = interp1(xv(ixrng), ts(ixrng), 0);             % ‘ts’ Values At Zero-Crossings
    plot(tzx{k1}, zeros(size(tzx{k1})), 'bp')               % Plot Zero Crossings
    hold off
end
xlabel('Time')

[Sam]

That helps so much! thanks :) So for the last problem, 3. How would I include the equals value in the equation set up you used for problem 1? Or is there a different format?

[Star Strider]

My pleasure!

For the particular solution, the function changes to:

Fp = @(t,x,w) [x(2); cos(t.*w)-x(1).^3.*0.2-x(1)-x(2).*0.2)];

and in the first loop, you change from iterating over ‘e’ to iterating over ‘w’, with the values for ‘w’ as stated in the problem. I will leave those changes to you, since they are straightforward. The second loop changes similarly to reflect the changes in the problem requirements.

[Sam]

Thanks! So , you just get everything to one side, that makes sense. Lastly, What exactly goes the green function represent? And if I wanted to make my graphs with just the graph of the function using the epsilon . Which part of the code would I take out to get rid of the green line.

[Star Strider]

My pleasure!

I’m not quite sure what you mean by ‘green function’. (I am using R2015a, and it uses the parula colormap as its default. The lines were blue and red when I plotted them.) In my code, I plotted both the function ‘x(:,1)’ and the first derivative, ‘x(:,2)’.

To plot only the function and not the derivative, just plot ‘x(:,1)’, or equivalently (in my latest code that also detects the first zero-crossing), the ‘xv’ variable. Since I saved them as cells, it will probably be easier to create and plot ‘xv’.

[saddam mollah]

Can you plot this problem in 3dimension as: t along x-axis, ev along y-axis and x1 along z-axis? Thank you in advance.