Help Needed for the following coding problem

So I have a problem at hand. I have few differential equations which are interdepent on one another. I have a written code with the following workflow :

1) I define intial conditions pretaining to the Vectors that I'll be using.

2) I run a for loop based on those on initial condtions for the differential equation

3) I apply Ode45 to the result to solve the differential equation and then send the values back to the starting point taking that as inital conditon and repeat this for say 40 steps.

I am facing the following errors -

1) My Ode45 solver doesn't work and I get 0 values as the updated final values.

2) If I go for a Multidimensional ode45 solver becasue I have like 8 D.E to be solved, it just doesn't work. I have attached a sample code and a sample image for reference.

Thanks for the help!

clc;
clear all;
close all;


N = 6; 


P = [1 2 1 2 3 1; 1 2 2 1 1 3];  



Distance = zeros(N,N);
for i = 1:N
    Distance(:,i) = sqrt(((P((2*i)-1)-P(1:2:(2*N-1))).^2)+((P(2*i)-P(2:2:(2*N))).^2)); %Distance between All Agents
end

Z = Distance < 3;%-eye(N)
A_undir = Z; %adjacency Matrix
D_undir = diag(sum(A_undir,1)); %Degree Matrix (Out-Degree)
L_undir = D_undir - A_undir; 



Raidus = 0.5;  
sigma = Raidus/2;
episilon = 0.5;
psi = 0.8;
epsilion_tilda = episilon + psi;
k1 = 1.5;
k2 = 0.15;
k3 = 1.5;
k4 = 0.0001;
kp = 6;
ki = 2;
k=5;
gamma = 1.5;


z1_dot = zeros(1,N);
z2_dot = ones(1,N);
w1_dot = zeros(1,N);
w2_dot = zeros(1,N);
z1 = zeros(1,N);
w1 = zeros(1,N);
z2 = zeros(1,N);
w2 = zeros(1,N);
Lambda2 = zeros(1,N);
Del_lambda_tilda = zeros(1,N);
V2_tilda = [1 1 1 1 1 1]; 
V2_tilda_dot = zeros(1,N);
alpha1 = V2_tilda;
alpha2 = (V2_tilda).^2;
A = A_undir;
U_c = zeros(2,N);
U_e = zeros(2,N);
P_dot = zeros(2,N);
Beta = Distance < 2;
Norm_Matrix = (Beta.^2)/(2*(sigma^2));

Po_x = zeros(N,N);
Po_y = zeros(N,N);


sumz1_1=z1.'-z1;
sumz1 = (sum(sumz1_1,2))';
sumw1_1=w1.'-w1;
sumw1 = (sum(sumw1_1,2))';
sumz2_2=z2.'-z2;
sumz2 = (sum(sumz2_2,2))';
sumw2_2=w2.'-w2;
sumw2 = (sum(sumw2_2,2))';

sumv2_tilda_1 = V2_tilda.'-V2_tilda;
sumv2_tilda = (sum(sumv2_tilda_1,2))';

sum_adjacency = (sum(A,2))';

for i = 1:N
            z1_dot(:,i) = gamma*(alpha1(:,i)-z1(:,i))-kp*sumz1(:,i)+ki*sumw1(:,i);
            w1_dot(:,i) = -ki*sumz1(:,i);
            z2_dot(:,i) = gamma*(alpha2(:,i)-z2(:,i))-kp*sumz2(:,i)+ki*sumw2(:,i);
            w2_dot(:,i) = -ki*sumz2(:,i);
            V2_tilda_dot(:,i) = -(k1*z1(:,i))-(k2*sum_adjacency(:,i)*sumv2_tilda(:,i))-(k3*(z2(:,i)-1)*V2_tilda(:,i))-(k4*abs(V2_tilda(:,i)).*V2_tilda(:,i));
            Lambda2(:,i) = k3*(1-z2(:,i))/k2;
            
            %Defining Position Matrix Based on P(i) - P(j)                       
            Po_x(:,i) =   (P((2*i)-1)-P(1:2:(2*N-1))) ;
            Po_y(:,i) = + (P(2*i)-P(2:2:(2*N)));  
            
            Del_lambda_tilda = -sum_adjacency(:,i)*sumv2_tilda(:,i)*([Po_x(i,:); Po_y(i,:)]).*(1/sigma)^2;  %Replace 1 with Position
            A = exp(-Norm_Matrix);               
            U_c(:,i) = ((csch(Lambda2(:,i)-epsilion_tilda)).^2)*Del_lambda_tilda(:,i);
            U_e(:,i) = [k*cos(2*pi/N+1)*i k*sin(2*pi/N+1)*i]';
            P_dot(:,i) = U_c(:,i) + U_e(:,i);     
end


t_initial = 0;
t_final = 100;
x_initial1 = zeros(1,N)';
x_initial2 = zeros(1,N)';
x_initial3 = zeros(1,N)';
x_initial4 = zeros(1,N)';
x_initial5 = ones(1,N)';
x_initial6 = P(end-1,:)';
x_initial7 = P(end,:)';

[t,z1_solver] = ode45(@(t,z1_solver) z1_dot',[t_initial t_final],x_initial1);
[t,w1_solver] = ode45(@(t,w1_solver) w1_dot',[t_initial t_final],x_initial2);
[t,z2_solver] = ode45(@(t,z2_solver) z2_dot',[t_initial t_final],x_initial3);
[t,w2_solver] = ode45(@(t,w2_solver) w2_dot',[t_initial t_final],x_initial4);
[t,V2_tilda_solver] = ode45(@(t,V2_tilda_solver) V2_tilda_dot',[t_initial t_final],x_initial5);
[t,P_dot_solver] = ode45(@(t,P_dot_solver) P_dot(end-1,:)',[t_initial t_final],x_initial6);
[t,P_dot1_solver] = ode45(@(t,P_dot1_solver) P_dot(end,:)',[t_initial t_final],x_initial7);


z1 = z1_solver(end,:);
w1 = w1_solver(end,:);
z2 = z2_solver(end,:);
w2 = w2_solver(end,:);
V2_tilda = V2_tilda_solver(end,:);
Position1 = P_dot_solver(end,:);
Position2 = P_dot1_solver(end,:);
P = [Position1 ; Position2];

Answer 1

[Walter Roberson]

z1_dot(:,i) = gamma*(1-z1(:,i))-kp*sumz1(:,i)+ki*sumw1(:,i);

Those create numeric values that will be effectively constants afterwards. They are not formulas in t or z1.

[t,z1_solver] = ode45(@(t,z1_solver) z1_dot',[t_initial t_final],x_initial1);

No matter what the t or z1_solver inputs, you always return the constants z1_dot'

Effectively you are integrating the z1' = constant over the time interval, using zeros (x_initial1) as the starting point. Integral of a constant is t(end) times the constant, minus t(1) times the constant, which is just a linear system.

[Pranjal Bhatia]

I have updated the code. The earlier one I posted was just a sample code. Am I still making the same mistake?

[Walter Roberson]

[t,z1_solver] = ode45(@(t,z1_solver) z1_dot',[t_initial t_final],x_initial1);

Yes. z1_dot is constant, so you are programming a linear system.

[Pranjal Bhatia]

Thanks for the help! Will try and solve this issue.

[Pranjal Bhatia]

You mentioned z1 being a constant. I am trying to fix that error, I have attached an image referring as to which equation I am trying to model in Matlab. I was assuming to be a first order ODE and using a loop to update the values. This is a dynamic average consensus estimator which z1 is trying to implement. I hope I am not asking too much at this point.

[Walter Roberson]

Unfortunately your posted code does not initialize a couple of variables such as w1, and your code seems to confuse ki and Ki. Because of these, I cannot tell what size your z1_dot will come out as. At the moment it looks like it will come out as a row vector, but your "in vector form" construction requires something that is 2 x 2 array.

Though if I understand correcly that your "in vector form" is hiding the fact that z and w are each vectors, then instead of something that is 2 x 2 array, then your z1_dot would need to be an array that is 200 x 200, to be matrix multiplied by the 200 x 1 array that will be the boundary conditions, and then a 200 x 1 vector added to it. The code would look like

[t,z1_solver] = ode45(@(t,zw) A*zw + B, [t_initial t_final], x_initial1);

where A is a 200 x 200 array, and B is a 200 x 1 array. zw will be a 200 x 1 vector of boundary conditions. (200 x 200) * (200 x 1) -> 200 x 1, add 200 x 1 to that is fine giving 200 x 1, so output would be 200 x 1 of derivative values, and that works fine in terms of sizes needed.

[Pranjal Bhatia]

I am sorry about that, while copying I made a mistake and that is why Ki and ki have been defined in two different ways, same thing regarding w1.

z1_dot is a row vector of 1XN here N has been described by me in the code. All are row vectors except P which defines the position.

would u suggest to describe z1_dot in a vector form rather than the way I defined it?

[Pranjal Bhatia]

Update: The code works, but what I want is that the final updated values are then used as intial values for the for loop. Can you help me with that?

[Walter Roberson]

Do not put the assignment of zeros to x_initial1 and related variables inside the loop if you want to use the updated values as the new initial points.

[Pranjal Bhatia]

Then how should I proceed? Change the inital values to the intial variables that were defined by me?

[Walter Roberson]

what I want is that the final updated values

Which variables hold those?

are then used as intial values for the for loop

Initial values for what variables?

What do you want to do about the fact that during the first run, there are no previous "final updated values" ?

Is the purpose of all of this that you are about to run (most of) the same code but with slightly different conditions, and you expect that the new outputs will be "close to" the results from the previous conditions, making it a good idea to start there?

Or is the idea that each time you loop through the same process, if you use the previous results, the new results will become more accurate, so you want to loop several times to improve accuracy ?

[Pranjal Bhatia]

Well for example, Look at z1_dot

1) First the for loop runs using the zero initial conditions described for z1_dot and z1

2) Then ode45 runs and gives me z1 from z1_dot.

3) Now I want the updated value of z1 to be used as an inital condition for z1_dot and then again continue this loop for a step number say 40.

The purpose is to find a varitaion , z1,w1,z2,w2 are dyanmic average PI consensus estimators and P_dot is the change in position.

[Walter Roberson]

Using the framework that you posted just for illustration purposes:

clc;
clear all;
close all;


N = 6; 


P = [1 2 1 2 3 1; 1 2 2 1 1 3];  



Distance = zeros(N,N);
for i = 1:N
    Distance(:,i) = sqrt(((P((2*i)-1)-P(1:2:(2*N-1))).^2)+((P(2*i)-P(2:2:(2*N))).^2)); %Distance between All Agents
end

Z = Distance < 3;%-eye(N)
A_undir = Z; %adjacency Matrix
D_undir = diag(sum(A_undir,1)); %Degree Matrix (Out-Degree)
L_undir = D_undir - A_undir; 



Raidus = 0.5;  
sigma = Raidus/2;
episilon = 0.5;
psi = 0.8;
epsilion_tilda = episilon + psi;
k1 = 1.5;
k2 = 0.15;
k3 = 1.5;
k4 = 0.0001;
kp = 6;
ki = 2;
k=5;
gamma = 1.5;

z1_dot = zeros(1,N);
z2_dot = ones(1,N);
w1_dot = zeros(1,N);
w2_dot = zeros(1,N);
z1 = zeros(1,N);
w1 = zeros(1,N);
z2 = zeros(1,N);
w2 = zeros(1,N);
Lambda2 = zeros(1,N);
Del_lambda_tilda = zeros(1,N);
V2_tilda = [1 1 1 1 1 1]; 
V2_tilda_dot = zeros(1,N);

A = A_undir;
U_c = zeros(2,N);
U_e = zeros(2,N);
P_dot = zeros(2,N);

t_initial = 0;
t_final = 100;
x_initial1 = zeros(1,N)';
x_initial2 = zeros(1,N)';
x_initial3 = zeros(1,N)';
x_initial4 = zeros(1,N)';
x_initial5 = ones(1,N)';


Po_x = zeros(N,N);
Po_y = zeros(N,N);

%initial values only above this point
    
for iterations = 1 : 10

    %below this point, values might be what was set after the end of
    %the previous iteration
    
    alpha1 = V2_tilda;
    alpha2 = (V2_tilda).^2;
    Beta = Distance < 2;
    Norm_Matrix = (Beta.^2)/(2*(sigma^2));
    
    x_initial6 = P(end-1,:)';
    x_initial7 = P(end,:)';
    
    sumz1_1=z1.'-z1;
    sumz1 = (sum(sumz1_1,2))';
    sumw1_1=w1.'-w1;
    sumw1 = (sum(sumw1_1,2))';
    sumz2_2=z2.'-z2;
    sumz2 = (sum(sumz2_2,2))';
    sumw2_2=w2.'-w2;
    sumw2 = (sum(sumw2_2,2))';
    
    sumv2_tilda_1 = V2_tilda.'-V2_tilda;
    sumv2_tilda = (sum(sumv2_tilda_1,2))';
    
    sum_adjacency = (sum(A,2))';
    
    for i = 1:N
                z1_dot(:,i) = gamma*(alpha1(:,i)-z1(:,i))-kp*sumz1(:,i)+ki*sumw1(:,i);
                w1_dot(:,i) = -ki*sumz1(:,i);
                z2_dot(:,i) = gamma*(alpha2(:,i)-z2(:,i))-kp*sumz2(:,i)+ki*sumw2(:,i);
                w2_dot(:,i) = -ki*sumz2(:,i);
                V2_tilda_dot(:,i) = -(k1*z1(:,i))-(k2*sum_adjacency(:,i)*sumv2_tilda(:,i))-(k3*(z2(:,i)-1)*V2_tilda(:,i))-(k4*abs(V2_tilda(:,i)).*V2_tilda(:,i));
                Lambda2(:,i) = k3*(1-z2(:,i))/k2;
                
                %Defining Position Matrix Based on P(i) - P(j)                       
                Po_x(:,i) =   (P((2*i)-1)-P(1:2:(2*N-1))) ;
                Po_y(:,i) = + (P(2*i)-P(2:2:(2*N)));  
                
                Del_lambda_tilda = -sum_adjacency(:,i)*sumv2_tilda(:,i)*([Po_x(i,:); Po_y(i,:)]).*(1/sigma)^2;  %Replace 1 with Position
                A = exp(-Norm_Matrix);               
                U_c(:,i) = ((csch(Lambda2(:,i)-epsilion_tilda)).^2)*Del_lambda_tilda(:,i);
                U_e(:,i) = [k*cos(2*pi/N+1)*i k*sin(2*pi/N+1)*i]';
                P_dot(:,i) = U_c(:,i) + U_e(:,i);     
    end
    
    [t,z1_solver] = ode45(@(t,z1_solver) z1_dot',[t_initial t_final],x_initial1);
    [t,w1_solver] = ode45(@(t,w1_solver) w1_dot',[t_initial t_final],x_initial2);
    [t,z2_solver] = ode45(@(t,z2_solver) z2_dot',[t_initial t_final],x_initial3);
    [t,w2_solver] = ode45(@(t,w2_solver) w2_dot',[t_initial t_final],x_initial4);
    [t,V2_tilda_solver] = ode45(@(t,V2_tilda_solver) V2_tilda_dot',[t_initial t_final],x_initial5);
    [t,P_dot_solver] = ode45(@(t,P_dot_solver) P_dot(end-1,:)',[t_initial t_final],x_initial6);
    [t,P_dot1_solver] = ode45(@(t,P_dot1_solver) P_dot(end,:)',[t_initial t_final],x_initial7);
    
    
    z1 = z1_solver(end,:);
    w1 = w1_solver(end,:);
    z2 = z2_solver(end,:);
    w2 = w2_solver(end,:);
    V2_tilda = V2_tilda_solver(end,:);
    Position1 = P_dot_solver(end,:);
    Position2 = P_dot1_solver(end,:);
    P = [Position1 ; Position2];
    
end

[Pranjal Bhatia]

Thanks a ton! You have been really helpful. Thank you :).

Just one last thing , Can i merge all of these ode45 together ? Like for reference : I know this won't work in this form as I am not calling ode45 correctly here. I was looking to solve this in a Multi Dimensional ode45

t0 = 0;
tf = 100;

Solver_Matrix = [ z1_dot w1_dot z2_dot w2_dot V2_tilda_dot P_dot(end-1,:) P_dot(end,:)]';
Initial_Conditions = [ zeros(1,N) zeros(1,N) zeros(1,N) zeros(1,N) ones(1,N) P(end-1,:) P(end,:)]'; 

tspan = [t0, tf];

[t, state_values] = ode45(Combined_Matrix,tspan,Initial_Conditions);

z1 = state_values(:,1);
w1 = state_values(:,2);
z2 = state_values(:,3);
w2 = state_values(:,4);
V2_Tilda = state_values(:,5);
Position1 = state_values(:,6);
Position2 = state_values(:,7);
P = [Position1 ; Position2];

[Walter Roberson]

If you have two ode with the same timespan, and the same options, and for which the same events would work (if events are configured), then you can merge them. You would concatenate the code and change the array offsets as needed.

Note: if you have the symbolic toolbox, then I would recommend that you have a look at configuring the odes as a system of symbolic equations, and that once you have done that, that you follow the work flow in the first example of odeFunction() in order to generate an anonymous function to pass to ode45()

[Pranjal Bhatia]

Sadly I don't have acess to Symbolic toolbox. Will go with the first method, the time span is the same only the initial conditions vary for each state equation. Thanks a lot for all your help :)