MATLAB ODE solver giving 0's everywhere

Question

My code is giving me zeros everywhere for my solution vectors but I don't know why. I've decomposed a coupled second order ODE into 4 1st order ODE's.

I have my function defined as xp.m

function zprime = f(t,z)
a = 1;
b = 1;
c = 1.5;

zprime = zeros(4,1);

zprime(1) = z(2);
zprime(2) = -a*z(1) + b*(z(3) - z(1));
zprime(3) = z(4);
zprime(4) = -c*(z(3) - z(1));
end

I run it in matlab with the following command:

[t,z] = ode45('xp',[1,100],[0 0 0 0]);

since my initial conditions are all 0. Is it that my initial conditions give the 0 solution or something else? When I change the ic's, the solutions change, as expected.

Thanks

Answer 1

Your actual question has more to do with mathematics than Matlab or programming. If you plug in zeros to your function f, you'll see immediately that it can give no other answer besides zero. You should look up equilibrium points or fixed points. Even if an equilibrium is unstable (imagine the top of a hill), a state exactly at that point with no disturbances (external inputs, numerical error) will remain balanced forever. If you're going to be working with differential equations, it's good to know how to find equilibrium points and then how to calculate the Jacobian matrix evaluated at these points in order to determine system properties. If you have further questions in this area, I'd suggest going to math.stackexchange.com.

Also, you're using an old scheme for calling your integration function. You can also pass in your parameters. Define your integration function as a sub-function in your main function or as a separate function file (same file name as function name – you'll want something other than f)

function zprime = f(t,z,a,b,c)
zprime(1,1) = z(2);
zprime(2,1) = -a*z(1) + b*(z(3) - z(1));
zprime(3,1) = z(4);
zprime(4,1) = -c*(z(3) - z(1));

Then, in your main function, call

a = 1;
b = 1;
c = 1.5;
[t,z] = ode45(@(t,z)f(t,z,a,b,c),[1 100],[0 0 0 0]);

Answer 2

For your particular case, with I.C.s z_0 = [0,0,0,0], the solution should be steady, with a value of z_out = [0,0,0,0]. Just take a look at your function, when you plug in z = z_0 and run it through your ODE solver.

zprime(1) = z(2);                       % --->  0
zprime(2) = -a*z(1) + b*(z(3) - z(1));  % ---> -a*(0) + b*(0)
zprime(3) = z(4);                       % ---> 0
zprime(4) = -c*(z(3) - z(1));           % ---> -c*(0-0) = 0

Keep in mind the general premise of a numerical solution. You take an initial condition, and feed it into the formula for the derivative. This tells you the slope of the function you are solving for. You use that slope to determine the value of the function at some future time (or nearby location, or similar) and then you feed that back into the derivative formula and start over again.

The only major difference between the different methods (Forward/backward Euler, RK4, ...) is the method that you use to determine the slope at the current location.