Errors when trying to plot a solution to a system of ODE

I try to plot on [1000, 5000] the solution of the system of ODEs

with the initial conditions , where and .

I used the function

function dzdt=odefunw1(t,z)
f=1/(t+1);
g=1+exp(-t);
h=diff(f);
dzdt=zeros(2,1);
dzdt(1)=z(2)-f*z(1);
dzdt(2)=(g+h)*z(1)-f*z(2);
end 

and the commands

tspan = [1000 5000];
z0 = [0.001 0.001];
[t,z] = ode45(@(t,z) odefunw1(t,z), tspan, z0);
plot(t,z(:,1),'r')

The following errors occured:

In an assignment  A(:) = B, the number of elements in A and B must be the same.

Error in odefunw1 (line 7)
dzdt(2)=(g+h)*z(1)-f*z(2);

Error in @(t,z)odefunw1(t,z)

Error in odearguments (line 90)
f0 = feval(ode,t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 115)
    odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

How could I fix it ? Many thanks in advance.

Answer 1

[Walter Roberson]

[t,z] = ode45(@(t,z) odefunw1(t,z), tspan, z0);

Okay, ode45 will invoke odefunw1

function dzdt=odefun1(t,z)

which will receive the parameters through the same names as in the outside (no complications about different variable names.)

ode45 always passes the first parameter, t, as a numeric scalar.

f=1/(t+1);

Using the / operator with a scalar on the right hand side is acceptable and will produce the same value as if you had used the ./ operator, so this line is okay

g=1+exp(-t);

Another scalar, not a problem

h=diff(f);

diff() applied to a numeric array is the numeric difference function that calculates the numeric difference between adjacent entries. You are passing it a numeric scalar. Because there are no adjacent entries to a scalar, the output of diff() applied to a numeric scalar is empty.

dzdt(2)=(g+h)*z(1)-f*z(2);

Because h is empty, the right hand side of that equation is empty.

When you look back at the mathematical formula, we see that your h corresponds to the formula term where . You should be taking the mathematic derivative of that, getting so inside the function you should calculate h as that formula.

[Cris19]

Thank you for the response. I fixed to

 dzdt=odefunw1(t,z)

I could take the derivative , but how could I use the derivative in the case when f is a piecewise function ? For example if one defines f as:

 if t >= 1
     f = 1/t;
   else
     f = -2*t^2+3*t;
   end;

what expression should be written in the line

   dzdt(2)=(g+h)*z(1)-f*z(2);

?

[Walter Roberson]

The rule about using piecewise functions with the ode*() series of functions is:

Don't

The ode*() series of functions use formulae that are only valid if all of the derivatives of the function that are used explicitly, plus two more derivatives (that is, Hessian is used), are continuous.

The sample definition of f you show is continuous at t = 1, but the second derivatives of it are not continuous at t = 1: f''(1) is -1 for the first form, but f''(1) = -4 for the second form. It is not valid to use ode*() with that piecewise definition of f in the case where t can cross 1.

In the case where the discontinuity is at a predictable time, it is advised to split the ode*() into two calls, one for the timespan up to the discontinuity, and one for the timespan starting from the discontuinuity.

In the case where the discontinuity is not at an easily predictable time (such as the case for a ball bounding), you would use event functions to trigger the end of ode*() processing, and then resume from the last time returned.

When you do the above discontinuity processing, if you are careful not to cross a discontinuity within any one invocation of ode*(), then you can program the different formulas in to be used conditionally. However in practice, doing so tends to mislead the programmer or readers of the code into thinking that if is generally valid inside the ode function, and then to faiing to take proper boundary conditions. Especially in the case where there is precisely one discontinuity at a predictable time, it can be clearer to write two different ode functions, one for each side of the discontinuity,

[t_left, z_left] = ode45(@left_side, [tspan(1), 1], z0);
z0_right = z_left(end,:);
[t_right, z_right] = ode45(@right_side, [1, tspan(2)], z0_right);
t = [t_left; t_right(2:end)];
z = [z_left; z_right(2:end,:)];
plot(t, z)

[Cris19]

Thank you so very much for the thorough answer!