CODEWITHSUNDEEP
matlabodeode45
Tina
Tina

Reputation: 23

Solving a system of ODEs using ODE45

I am trying to learn how to use MATLAB to solve a system of differential equations (Lorenz equations) and plot each solution as a function of t

X’ = −σx + σy  
Y’ = ρx − y − xz
Z’ = −βz + xy

where σ = 10, β = 8/3, and ρ = 28, as well as x(0) = −8, y(0) = 8, and z(0) = 27.

Here is the code that I am using:

function xprime = example(t,x)

sig    = 10;
beta   = 8/3;
rho    = 28;
xprime = [-sig*x(1) + sig*x(2); 
           rho*x(1) - x(2) - x(1)*x(3); 
          -beta*x(3) + x(1)*x(2)];

x0    = [-8 8 27]; 
tspan = [0 20]; 
[t,x] = ode45(@example, tspan, x0);

figure 
plot(t,x(:,1)), hold on
plot(t,x(:,2)), hold on
plot(t,x(:,3)), hold off

However, this generates an error, how do I fix this? I am not sure what input arguments are missing or where I am going wrong. I appreciate any help, thanks.

Not enough input arguments.

Error in example (line 9) xprime=[-sigx(1) + sigx(2); rho*x(1) - x(2) - x(1)x(3); -betax(3) +
x(1)*x(2)];

Upvotes: 2

Views: 2097

Answers (1)

Rody Oldenhuis
Rody Oldenhuis

Reputation: 38032

That is actually a pretty good first try!

The problem is that when you press the Run button (or press F5), you're calling the function example with no arguments; which is what MATLAB is complaining about.

A second problem is that, even if you were to be able to run the function like this, ode45 would call the function example, which would call ode45, which would call example, which would call ode45 and so on, until the recursion limit is reached.

The solution to both is to split it up in two functions (these may be written into the same M-file): 

% top-level function; no arguments
function caller() 

    x0    = [-8 8 27]; 
    tspan = [0 20]; 
    [t,x] = ode45(@example, tspan, x0);

    figure 
    plot(t,x(:,1)), hold on
    plot(t,x(:,2)), hold on
    plot(t,x(:,3)), hold off

end

% The derivative 
function xprime = example(t,x)

    sig    = 10;
    beta   = 8/3;
    rho    = 28;
    xprime = [-sig*x(1) + sig*x(2); 
               rho*x(1) - x(2) - x(1)*x(3); 
              -beta*x(3) + x(1)*x(2)];              
end

Upvotes: 3

Related Questions