Solving an ODE with changing boundary conditions

I have written a function to solve an ODE used to simulate a swarm of robots surrounding a moving point. What made sense to me was to have the initial and boundary conditions updated each second, and use a for loop to do this. The boundary conditions represent some of the robots which have knowledge of where the target curve/point is because of better sensing capabilities.

Parameters set:

%% Model Parameters

a= .1; % connection strength

r = 1; % radius of circle

L=.2;

f=@(t,z) 0.01; % simplified model dynamics

k=1.5; % boundary feedback gain

gamma=@(theta, t) [r*cos((2*pi)*theta);r*sin((2*pi)*theta)] + [1*t; 1*t]; % target curve

% tauM=.89; % time-delay bound

alpha=.1;

%% Simulation parameters

tf=20; % simulation time

N=20; % Number of agents

The function returns the error in x and y relative to the target curve (y0), over timespan tmesh0.

function [tmesh0,y0]=Numeric_soln_ODEs_const(a,f,k,gamma,tf,N)

% ODE modelling of MAS with a delay in the leader

disp('Solving ODEs with constant piecewise global controller...')

tic

deltaf=@(t,gamma,y) f(t,gamma+y)-f(t,gamma);

%% First iteration (t=0)

leader=ode45(@(t,y) -k*y+deltaf(t,gamma(0,0),y),[0,1],-gamma(0,0));

% Followers

h=1/(N+1);

T=-a/h*(diag(ones(1,N))-diag(ones(1,N-1),-1)); % Laplacian matrix

I_n=eye(size(gamma(0,0),1)); % n-dimensional identity matrix

gammaVal= gamma((1:N)*h,0); % Target position for the followers

disp("gammaVal size = " + num2str(size(gammaVal)));

rhs=@(t,y) kron(T,I_n)*y+deltaf(t,gammaVal(:),y)+kron([a/h;zeros(N-1,1)],I_n)*deval(leader,t);

[tmesh0,y0]=ode45(rhs,[0,1],-gammaVal(:));

% Combining data

y0=[deval(leader,tmesh0)',y0];

[a,b]=size(y0);

disp("rows="+num2str(a)+" cols="+num2str(b));

x1=gammaVal(1,:)+y0(:,3:2:end); % x values of agents

x2=gammaVal(2,:)+y0(:,4:2:end); % y values of agents

next_init = [x1(end,:);x2(end,:)]; % next iteration is the last spatial coordinates of the agents

disp("next_init size = " + num2str(size(next_init)));

%% Loop

for i = 1:1:tf-1

% Leader (Boundary conditions)

leader=ode45(@(t,y) -k*y+deltaf(t,gamma(0,i),y),[i,i+1],gamma(0,i));

% Followers

gammaVal= gamma((1:N)*h,i); % Target position for the followers

rhs=@(t,y) kron(T,I_n)*y+deltaf(t,gammaVal(:),y)+kron([a/h;zeros(N-1,1)],I_n)*deval(leader,t);

[tmesh,y]=ode45(rhs,[i,i+1],next_init);

% Combining data

y=[deval(leader,tmesh)',y];

y0 = [y0; y];

tmesh0 = [tmesh0; tmesh];

[a,b]=size(y0);

disp("rows="+num2str(a)+"cols="+num2str(b));

x1=gammaVal(1,:)+y0(:,3:2:end); % x values of agents. Have removed a value so matrix sizes match

x2=gammaVal(2,:)+y0(:,4:2:end); % y values of agents. Have removed a value so matrix sizes match

next_init = [x1(end,:);x2(end,:)]; % next iteration is the last spatial coordinates of the agents

end

disp("ODEs with global constant piecewise controller are solved in " + num2str(toc) + " sec")

Then to find the L2-Norm, I use this code in my main script:

figure('name','L2-norms of the error') ;

plot(Tode,sum(Yode.^2,2)/(N+1));

When I plot this, the error is of order 9, which is large enough for me to realise I have made a mistake when formulating. Anyone know where this could be, or just a simpler way to do it as I have a feeling I am overcomplicating?

[Torsten]

Do you think anybody will be able to understand what you are doing in your code without supplying a proper mathematical description of the problem ?

[Elizabeth Violet]

My main question is is there a way to include changing boundary conditions for an ODE in MATLAB, preferably using ode45? I'm not expecting you to understand everything that's going on in the code, its just to provide a vague background of how I'm currently doing it.

Currently, I'm using a for loop to update my initial conditions to be the last values of the previous iteration, and my boundary conditions to be the target curve at that moment. My target curve is set up as a circle moving away from the origin by 1m/s horizontally and vertically. I'm pretty sure this method is giving me te wrong answers, but I believe my method is pretty shoddy, even if it were giving me a more believable answer.

Having known boundary conditions changing at a known rate seems like something MATLAB might provide resources for, but I've not found them - would anyone know a function or way of passing a time variable that could be useful?

[Torsten]

% Leader (Boundary conditions)
leader=ode45(@(t,y) -k*y+deltaf(t,gamma(0,i),y),[i,i+1],gamma(0,i));    
% Followers
gammaVal= gamma((1:N)*h,i); % Target position for the followers 
rhs=@(t,y) kron(T,I_n)*y+deltaf(t,gammaVal(:),y)+kron([a/h;zeros(N-1,1)],I_n)*deval(leader,t); 
[tmesh,y]=ode45(rhs,[i,i+1],next_init); 

You integrate the leader and follower equations subsequently in steps of 1. Wouldn't it be possible to integrate both in one call to ode45 over the complete timespan [1 tf] ? Would this solve your "boundary conditions problem" (which I still don't understand) ?

[Elizabeth Violet]

It would be, and this is exactly what I do for an earlier iteration of the code which has a static target curve.

My issue is that my target curve is now moving, so the boundary conditions are changing constantly. I've coded it in a for loop which is definitely not the best way. I initiially tried passing a time variable to the condition by:

leader=ode45(@(t,y) -k*y+deltaf(t,gamma(0,t),y),[0 tf],gamma(0,t));

% Followers

gammaVal= gamma((1:N)*h,t); % Target position for the followers

rhs=@(t,y) kron(T,I_n)*y+deltaf(t,gammaVal(:),y)+kron([a/h;zeros(N-1,1)],I_n)*deval(leader,t);

[tmesh,y]=ode45(rhs,[0, tf],gamma(0,0));

But this throws an error to do with t being undefined, despite it being defined as a variable in each case it is used.

My formulation of the problem is relatively difficult to understand, but I'll go through a simpler example which has the same boundary condition problem:

Imagine a rod being heated on one end. The temperature at the very end of this rod is always known but also always changing, and there is no upper bound to this value (ie as the time goes to infinity, so too does the temperature at the end of this rod). The rate of change is also known and is a constant. How should I formulate the boundary conditions in MATLAB to take this into account?

[Torsten]

How should I formulate the boundary conditions in MATLAB to take this into account?

Assuming the boundary condition changes at the end of the rod and you formulate the temperature distribution of the rod as a PDE

rho*cp*dT/dt = d/dx(lambda*dT/dx)

and discretize this PDE in xstart = x1 < x2 < ... < xn = xend to get the temperatures T1, T2,..., Tn, the differential equation in xend for T would be

dTn/dt = known rate