Modelling a system of differential equations with recurrences in matlab

Trying to model a system in the form

a*u_n(t)'' + b*u_n(t) = k( v_n+1(t) + v_n-1(t) - 2u_n(t) )

c*v_n(t)'' + d*v_n(t) = k( u_n+1(t) + u_n-1(t) - 2v_n(t) )

a,b,c,d,k are all constants

Pretty sure this can only be done numerically

[Torsten]

This system can only be solved for a finite start and end value for n.

So you need "boundary conditions" that prescribe u_0 and v_0 ("start functions") as well as u_N and v_N ("end functions").

Further initial conditions at t = 0 are required for u_1,...,u_N-1 and v_1,...,v_N-1 as well as u_1',...,u_N-1' and v_1',...,v_N-1'.

Could you supply both boundary and initial conditions ?

[David Goodmanson]

Hi Austin, does the last term on the first line equal -2u_n(t) (certainly possible) or should it be -2v_n(t) ? Similarly for the last term on the second line.

[Austin]

The last terms of each equation are correct as is

For initial conditions u_i(0) = v_i(0) = 0 i=1..n

There are equations for what would be the boundary conditions

a*u_0(t)'' + b*u_0(t)=k( v_1(t) - u_0(t) )

c*v_0(t)'' + d*v_0(t)=k( u_1(t) - v_0(t) )

and

a*u_n(t)'' + b*u_n(t)=k( v_n-1(t) - u_n(t) )

c*v_n(t)'' + d*v_n(t)=k( u_n-1(t) - v_n(t) )

If its not possible from these then setting u_0, v_0, u_n, v_n to be any arbitrary constant/function might work well enough

[Torsten]

Since the boundary conditions are defined by second-order differential equations for u_0, v_0, u_n and v_n, we need u_i(0), u_i'(0), v_i(0), v_i'(0) for i = 0,...,n.

You said we may assume u_i(0) = v_i(0) = 0 for i=1,...,n-1. So u_0(0),v_0(0),u_n(0),v_n(0) and all derivatives u_i'(0) and v_i'(0) at t = 0 ( (i = 0,...,n) are to be added to the problem description to make the system solvable.

[Austin]

Of course, the i should be i=0..n not 1..n

Then for the derivatives u_i(0)' = 0.5, v_i(0)' = 0.3

Answer 1

[Torsten]

N = 20;

k = 1;
a = 1;
b = 1;
c = 1;
d = 1;

u0    = zeros(N,1);
udot0 = 0.5*ones(N,1);
v0    = zeros(N,1);
vdot0 = 0.3*ones(N,1);
y0 = [u0;udot0;v0;vdot0];

tspan = linspace(0,10,20);

[T,Y] = ode15s(@(t,y)fun(t,y,N,k,a,b,c,d),tspan,y0);

u    = Y(:,1:N);
udot = Y(:,N+1:2*N);
v    = Y(:,2*N+1:3*N);
vdot = Y(:,3*N+1:4*N);

function dydt = fun(t,y,N,k,a,b,c,d)
  u    = y(1:N);
  udot = y(N+1:2*N);
  v    = y(2*N+1:3*N);
  vdot = y(3*N+1:4*N);
  
  dudt   = zeros(N,1);
  d2udt2 = zeros(N,1);
  dvdt   = zeros(N,1);
  d2vdt2 = zeros(N,1);

  dudt = udot;  
  dvdt = vdot;

  %a*u_1(t)'' + b*u_1(t) = k( v_2(t) - u_1(t) )
  %a*u_i(t)'' + b*u_i(t) = k( v_i+1(t) + v_i-1(t) - 2u_i(t) ) (2 <= i <= N-1)
  %a*u_N(t)'' + b*u_N(t) = k( v_N-1(t) - u_N(t) )
  d2udt2(1)     = (-b*u(1)+k*(v(2)-u(1)))/a;
  d2udt2(2:N-1) = (-b*u(2:N-1)+k*(v(3:N)+v(1:N-2)-2*u(2:N-1)))/a;
  d2udt2(N)     = (-b*u(N)+k*(v(N-1)-u(N)))/a;
  
  %c*v_1(t)'' + d*v_1(t) = k( u_2(t) - v_1(t) )
  %c*v_i(t)'' + d*v_i(t) = k( u_i+1(t) + u_i-1(t) - 2v_i(t) ) (2 <= i <= N-1)  
  %c*v_N(t)'' + d*v_N(t) = k( u_N-1(t) - v_N(t) )  
  d2vdt2(1)     = (-d*v(1)+k*(u(2)-v(1)))/c;
  d2vdt2(2:N-1) = (-d*v(2:N-1)+k*(u(3:N)+u(1:N-2)-2*v(2:N-1)))/c;
  d2vdt2(N)     = (-d*v(N)+k*(u(N-1)-v(N)))/c;
  
  dydt = [dudt;d2udt2;dvdt;d2vdt2];
end

[Austin]

Thanks this was my first time encountering a system of equations like this so this helped a lot

Answer 2

[John D'Errico]

This is known as a delay differential equation. You will find any solvers for them starting with the letters dde.

help dde23

 dde23 - Solve delay differential equations (DDEs) with constant delays
    This MATLAB function, where tspan = [t0 tf], integrates the system of
    delay differential equations y′(t)=f(t,y(t),y(t−τ1),...,y(t−τk)) over
    the interval specified by tspan, where τ1, ..., τk are constant,
    positive delays specified by delays.

    Syntax
      sol = dde23(ddefun,delays,history,tspan)
      sol = dde23(ddefun,delays,history,tspan,options)

    Input Arguments
      ddefun - System of delay differential equations to solve
        function handle
      delays - Time delays
        positive vector
      history - Solution history
        function handle | vector | structure
      tspan - Interval of integration
        two-element vector
      options - Integrator options
        structure array

    Output Arguments
      sol - Solution for evaluation
        structure array

    Examples
      openExample('matlab/DDE23ConstantDelaysExample')
      openExample('matlab/LocateZeroCrossingsOfDDEExample')

    See also ddesd, ddensd, ddeget, ddeset, deval

    Introduced in MATLAB before R2006a
    Documentation for dde23
       doc dde23

You will convert the second order DDEs each into a pair of first order DDEs using the standard trick, so you will have a system of 4 DDEs. Standard trick:

If you have a second order equation of the form:

y''(x) = stuff

you convert it into a pair of first order equations by creating a new unknown function, I'll call it z, where z is just the currently unknown first derivative of y.

y'(x) = z(x)

z'(x) = stuff

The same will apply in your case, even with a DDE.

[Austin]

Sorry my system of equations were written out poorly each u_i, v_i i=1..n+1 is a funtion of t and has been updated in question. Not sure if dde would still be used in this case