How to solve a system of ODEs and plot the result

Hello, I am fairly new to Matlab. I am trying to use matlab to solve an ODE with known boundary conditions and plot the results against those I've obtained analytically. The ODE in question is

d^2y/dt^2 + 400y = 0

For t 0-5 seconds and initial displacement y(0)=0.1 and velocity dy/dt=0. Since it has a second derivative, my understanding from reading so far is that I need to rewrite it as a sytem of equations with only first order derivatives. So subsituting y(1)=y(t) and y(2)=dydt, I tried using this code that I found,

function dydt = odefun(t,y)

dydt = zeros(2,1);

dydt(1) = 0*y(1)+y(2);

dydt(2) = -400*y(1)+0*y(2);

But it says that I don't have enough input arguements. I added in the 0* terms in hopes of correcting that, but no luck. I'm thinking maybe its asking for a t domain, but I tried a couple ways of defining that without sucess. I also played around with a couple alternative codes but no dice there either. I'm thinking if I get this ODE function to work, I insert it into ODE45 solver and plot the results. Thanks in advance for any help.

Answer 1

[James Tursa]

This is where you needed to show us the complete code and the complete error message. It is probably complaining about your initial conditions input, but without seeing your complete code I can't be sure. So something like this might be the fix:

[T,Y] = ode45(@odefun,[0 5],[0.1;0]); % Two-element initial conditions

Adding those 0* terms in your derivative function has no effect and can be removed.

[Thomas Ward]

I think that worked, I didn't realize it needed the ode45 part to be considered fully defined, I thought I could troubleshoot this part of the code and then add the ODE solver at the end. The entirety of my code is now:

[T,Y] = ode45(@odefun,[0 3],[0.1;0]);
t=linspace(0,3);
undampedx = 0.1 * cos(10*t);

figure(1)
plot(t,undampedx,T,Y)
function dydt = odefun(t,y)
dydt = zeros(2,1);
dydt(1) = 0*y(1)+y(2);
dydt(2) = -400*y(1)+0*y(2);
end

[Tanmayi Pagadala]

what does the "dydt = zeros(2,1);" mean?

[Walter Roberson]

@Tanmayi Pagadala

Suppose that line were omitted -- suppose the code was

function dydt = odefun(t,y)
dydt(1) = 0*y(1)+y(2);
dydt(2) = -400*y(1)+0*y(2);
end

then that first line dydt(1) = 0*y(1)+y(2); would create dydt as a scalar variable and assign the value on the right hand side of the equation. After that line, dydt will exist in the workspace (it does not exist before that assignment.)

Then you would encounter the line dydt(2) = -400*y(1)+0*y(2); which asks to extend the scalar variable dydt to a second location and assign the second location the value on the right hand side.

But... when it extends that scalar to a second memory location, where is that second memory location relative to the first memory location? Does it become the second row of a 2 x 1 column vector? Does it become the second column of a 1 x 2 row vector? Somewhere else?

The answer is that in MATLAB if you have a scalar and you ask to assign additional elements to it, and use a single index to do that -- not dydt(2,1) or dydt(1,2) but just dydt(2) -- then MATLAB always creates a row -- so dydt would become 1 x 2.

However... it happens that ode45() and similar functions require your ode function to return a column vector -- 2 x 1 rather than 1 x 2 in this case.

You can handle that in three ways:

1. After you assign to all elements of dydt (getting a row vector), you can assign dydt = dydt.' to transpose the 1 x 2 into a 2 x 1; OR
2. After you assign to all elements of dydt (getting a row vector), you can reshape() it into a column vector using either dydt = reshape(dydt, [], 1) or dydt = dydt(:) ; OR
3. Before you assign anything to dydt(1) and dydt(2) you can initialize dydt to the size and shape you want it to end up with using dydt = zeros(2,1) . When you have an existing 2 x 1 array and you ask to assign to dydt(2) then MATLAB looks and sees that the array already exists and is large enough and puts the value into the column; you do not need to transpose or reshape or (:) afterwards.