Using ODE, several errors with a spring system and function

1 visualización (últimos 30 días)
Julia Russo
Julia Russo el 12 de Jul. de 2016
Respondida: Mischa Kim el 13 de Jul. de 2016
I am attempting to set up a function file for the following differential equation: (u''+ u + epsilon*u^3 = 0) with initial conditions: u(0) = 0, u'(0) = 1.
My variables are t, epsilon, and u, so my code for my function is:
function [ du ] = fu( t,e,u )
du= zeros(2,1);
du(1)= u(2);
du(2)= -(u(1))-e*(u(1))^3;
end
saved as 'F1'
when I run:
e=[0.0:.2:1.0]
t=[1 20]
ic=[0 1]
[T, U]=ode45('F1',t,e,ic)
I get a huge series of errors starting with: Attempted to access u(2); index out of bounds because numel(u)=0.
Error in F1 (line 3)
du(1)= u(2);
Please help me! I have no clue why this is not working and I'm extremely stuck.
  2 comentarios
James Tursa
James Tursa el 12 de Jul. de 2016
What is the purpose of e being a vector? Are you trying to solve the ODE multiple times for multiple values of epsilon? Do you need help with how to pass that epsilon value into your F1 routine?
Robert Fennis
Robert Fennis el 12 de Jul. de 2016
Have you saved all of your latest files? Also make sure the function fu has the same name as your file. I don't know if this causes the problem but you are not supposed to mix those up. Also, I see that they use a different syntax in the examples:
[t,y] = ode45(@vdp1,[0 20],[2; 0]);

Iniciar sesión para comentar.

Respuestas (1)

Mischa Kim
Mischa Kim el 13 de Jul. de 2016
Julia, use
function myODE()
e = 0.0:.2:1.0; % no brackets, the colons are used to define the vector
t = [1 20];
ic = [0; 1]; % column (not row) vector
for ii = 1:numel(e) % I guess you want to solve for different e?
[T, U] = ode45(@fu,t,ic,[],e(ii));
hold on
plot(T,U)
end
end
function [du] = fu(~,u,e)
du = [u(2); -(u(1))-e*(u(1))^3]; % like above, column not row vector
end
Put the code above into one .m file, name it myODE.m and run.

Categorías

Más información sobre Ordinary Differential Equations en Help Center y File Exchange.

Community Treasure Hunt

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

Start Hunting!

Translated by