Help Needed for the following coding problem

1 vue (au cours des 30 derniers jours)
Pranjal Bhatia
Pranjal Bhatia le 3 Juil 2020
Commenté : Pranjal Bhatia le 4 Juil 2020
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];

Réponses (1)

Walter Roberson
Walter Roberson le 3 Juil 2020
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.
  15 commentaires
Walter Roberson
Walter Roberson le 4 Juil 2020
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
Pranjal Bhatia le 4 Juil 2020
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 :)

Connectez-vous pour commenter.

Tags

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by