How to solve a system of ODEs and plot the result

19 visualizaciones (últimos 30 días)
Thomas Ward
Thomas Ward el 10 de Feb. de 2020
Comentada: Walter Roberson el 23 de Mzo. de 2022
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.

Respuesta aceptada

James Tursa
James Tursa el 10 de Feb. de 2020
Editada: James Tursa el 10 de Feb. de 2020
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.
  3 comentarios
Tanmayi Pagadala
Tanmayi Pagadala el 23 de Mzo. de 2022
what does the "dydt = zeros(2,1);" mean?
Walter Roberson
Walter Roberson el 23 de Mzo. de 2022
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.

Iniciar sesión para comentar.

Más respuestas (0)

Categorías

Más información sobre Programming en Help Center y File Exchange.

Etiquetas

Productos


Versión

R2018b

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by