solve ode with ode45

Hello,

So I have three different parameters for dr ( 0.2, 1, 2) and wanted to know the proper way to pass it through the ode equation without an error.

with initial condition x(0) = 2 and x'(0) = 0 and over the time interval 0<t<5

[t,y] = ode45(@odepractice,t,[x0 v0]);
xsim = y(:,2); 
vsim = y(:,n); 

function dy = odepractice(t,y)
dr = 0.2 1 2
wp = 5
dy = zeros(n,2);

dy(1) = y(2);
dy(2) = -2*dr*wp*y(2) - wp^2*y(1);
end

Error using odearguments

ODEEXAMPLE returns a vector of length 5, but the length of initial conditions vector is

2. The vector returned by ODEpractice and the initial conditions vector must have the

same number of elements.

Error in ode45 (line 107)

odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in Problem2 (line 87)

[t,y] = ode45(@odepractice,t,[x0 v0,]);

Error in Problem2 (line 87)

[t,y] = ode45(@odepractice,t,[x0 v0,]);

[Rena Berman]

(Answers Dev) Restored edit

Answer 1

[Sam Chak]

Hi @Lan

If I'm not mistaken, 'dr' stands for 'damping ratio.' You probably want to observe how the response dampens over time. You should try out the approaches suggested by @Walter Roberson. For a simple problem, I use the for-loop approach.

dr    = [0.2 1.0 2.0];  % different values of damping ratio (zeta)
x0    = 2;
v0    = 0;
y0    = [x0 v0];
t     = [0   6];

for j = 1:numel(dr)     % loop 3 times because there are 3 zeta values
    [t, y] = ode45(@(t, y) odepractice(t, y, dr(j)), t, y0);
    xsim   = y(:,1); 
    vsim   = y(:,2);
    plot(t, xsim), hold on
end

hold off, grid on
xlabel('t'), ylabel('x(t)')
legend('\zeta = 0.2', '\zeta = 1.0', '\zeta = 2.0', 'fontsize', 12)

function dy = odepractice(t, y, dr)
    wp    =   5;
    dy    =   zeros(2, 1);
    dy(1) =   y(2);
    dy(2) = - 2*dr*wp*y(2) - wp^2*y(1);
end

[Lan]

So when I tried to overlay the xsim with my other anlaytical solution, I got an error saying the vectors must be the same length. The length of the analytical solution is 1000, but the length for xsim is 157. How can I make the length of xsim to be the same length as the analytical (1000)?

[Torsten]

hold on
plot(tsim,xsim)
plot(tanalytic,xanalytic)
hold off

[Lan]

I meant to say changing the vector length of (t/xsim) to a 1000. Iis there a way to write a code to fix that criteria? thanks

[Sam Chak]

Hi @Lan

You can try out @Torsten's suggestion and see what happens. Due to the relatively high accuracy of the numerical Runge–Kutta solver, ode45(), it is highly likely that you won't observe any significant differences between the numerical solution and the analytical solution in the same graph. Therefore, I suggest calculating the error and plotting it as shown below:

t         = linspace(0, 3, 3001);
y0        = [2 0];
[tsim, y] = ode45(@odepractice, t, y0);
xsim      = y(:,1); 

% Plot Comparison
subplot(2,1,1)
% plot numerical solution
plot(tsim, xsim), hold on

% compute analytical solution based on generated tsim
xanalytic = 2*exp(- 5*tsim).*(1 + 5*tsim);  % <-- equation depends on tsim
plot(tsim, xanalytic), hold off, grid on

xlabel('t'), ylabel('x(t)')
title('Comparison between x_{sim} and x_{anal}')
legend('x_{sim}', 'x_{anal}', 'fontsize', 12)

% Plot Error
subplot(2,1,2)
error     = xanalytic - xsim;
plot(tsim, error), grid on
xlabel('t'), ylabel('e(t)')
title('Error between x_{sim} and x_{anal}')

function dy = odepractice(t, y)
    wp    =   5;
    dr    =   1;
    dy    =   zeros(2, 1);
    dy(1) =   y(2);
    dy(2) = - 2*dr*wp*y(2) - wp^2*y(1);
end

Answer 2

[Walter Roberson]

dr = 0.2 1 2

That is not valid MATLAB syntax.

n is not defined.

Your dr looks like it is trying to be a vector of length 3, so when you use it in defining dy(2) it looks like you would get an error that the right side is a different length than the left side.

It is not clear what your intent is for having 3 different values for dr.

[Walter Roberson]

http://www.mathworks.com/help/matlab/math/parameterizing-functions.html

[Lan]

The intent is to use that values to plot anlaytical and numerical graph.

n is suppose to be a function of t

[Walter Roberson]

If you want to run an ode with the same equation, but several different values of a constant, there are three possible ways to do it:

1. You can construct several different versions of the ode function, one version for each value of the constant, and run each version individually;
2. You can construct a single version of the ode function that accepts the "current" value of the constant as a parameter, and loop over the permitted values of the parameter passing each different value in turn in an ode45() call using "parameterization" of the anonymous function (instructions linked to above); or
3. You can construct a single version of the function that (number of basic states) times (number of constant values) states, and run them all simultaneously in a single ode45() call.

For the last possibility, you can use vectorized calculations, something like

dy(1:2:end) = y(2:2:end);
dy(2:2:end) = -2*dr(:).*wp.*y(2:2:end) - wp.^2.*y(1:2:end);

Answer 3

[Torsten]

Comment out

t = [0 5]

before calling ode45 in the loop:

%t = [0 5];

[Lan]

I'm still getting this error:

Error using plot

Vectors must be the same length.

Error in Problem2 (line 103)

plot(t,x1/cm2m,'k',t,xtrig1/cm2m,'b',t,xsim/cm2m,'r:'); hold on

[Walter Roberson]

for j = 1:numel(df)
    [t, y] = ode45(@(t,y) odeexample(t,y,df(j)), t,y0);
    xsim = y(:,1); 
    vsim = y(:,2); 
end

Each time through the loop, you are overwriting all of t and all of y, and all of xsim and vsim . The only point in running that loop (instead of doing just the final iteration) is that if the ode45 call ends early (because it fails to integrate) then the t it returns will cover only the times up to when it stopped integrating, and then that shorter t will be passed to the next iteration's ode45() call -- so t can get shorter in the loop but not longer. It should not be relied upon that after that loop that might be shortening t that the final t will be the same size as the original t

Because the final xsim and vsim are being completely overwritten, the only purpose of doing those loops is to find the shortest time interval for which ode45() can integrate all of the different df values.

[Torsten]

@Lan

I'm still getting this error:

Error using plot

Vectors must be the same length.

Then you didn't make the change to your code I suggested or you use a different code (look above to see that it works). Or you chose a different value for cm2m.