I have given matrix and function to solve ode45 .here 'x' and 'f' is a column 20*1 vector. xprime is a derivative of x.please help me the code to complete .

##matrix M is here,

clc;

clear all;

n = 20;

C= zeros(n, n);

C(1,n)=1;

C(1,n-1)=1;

C(2,n)=1;

A = zeros(n, n);

for i = 1:n

for j = 1:n

if j==i+1 || j == i + 2

A(i,j) = 1;

else

A(i, j) = 0;

end

end

end

B1=A+C;

B=B1+B1';

% Find non-zero and zero indices in the upper triangular part of the matrix

nonZeroIndices = []; % Initialize an array to store non-zero indices

zeroIndices = []; % Initialize an array to store zero indices

for i = 1:n

for j = i+1:n

if B1(i, j) ~= 0

nonZeroIndices = [nonZeroIndices; i, j];

else

zeroIndices = [zeroIndices; i, j];

end

end

end

% Shuffle the non-zero and zero indices randomly

shuffledNonZeroIndices = nonZeroIndices(randperm(size(nonZeroIndices, 1)), :);

shuffledZeroIndices = zeroIndices(randperm(size(zeroIndices, 1)), :);

% Choose the first two non-zero and zero indices from the shuffled arrays

numChosen = 2;

chosenNonZeroIndices = shuffledNonZeroIndices(1:numChosen, :);

chosenZeroIndices = shuffledZeroIndices(1:numChosen, :);

% % Define values to replace the chosen positions

% replacementValueNonZero = 999; % Replace chosen non-zero positions with this value

% replacementValueZero = -999; % Replace chosen zero positions with this value

% Replace the chosen non-zero positions with the replacement value

for i = 1:2

B1(chosenNonZeroIndices(i, 1), chosenNonZeroIndices(i, 2)) = 0;

end

% Replace the chosen zero positions with the replacement value

for i = 1:2

B1(chosenZeroIndices(i, 1), chosenZeroIndices(i, 2)) = 1;

end

% Display the modified matrix B

disp("Modified Matrix B1:");

disp(A);

M=B1+B1';

disp(M)

%G=graph(M)

%total=nnz(M);

%plot(G)

%%

function

function xprime = fun(t, x)

theta = 0.1;dh = 0.2;phi = 0.3;ita = 0.15;

nn=20;

i = 1; f = [];

while i< nn

f = [f; x(i)*(1 - theta*x(i)) - (x(i)*x(i+1))/(1 + x(i));(phi*x(i)*x(i+1))/(1 + x(i)) - ita*x(i+1)];

i = i + 1;

end

xprime = f+dh*K*x;

end

%%

x0 = rand(20,1);

tspan=[0 1000];

[t,x] = ode45(funnn,tspan, x0);

[Star Strider]

The function needs to have ‘K’ passed to it, however I cannot determine what ‘K’ is here. The function also needs to be at the end of the code, after the ode45 call.

n = 20;
C= zeros(n, n);
C(1,n)=1;
C(1,n-1)=1;
C(2,n)=1;
A = zeros(n, n);
for i = 1:n
    for j = 1:n
        if j==i+1 || j == i + 2
            A(i,j) = 1;
        else
            A(i, j) = 0;
        end
    end
end
B1=A+C;
B=B1+B1';
% Find non-zero and zero indices in the upper triangular part of the matrix
nonZeroIndices = [];  % Initialize an array to store non-zero indices
zeroIndices = [];     % Initialize an array to store zero indices
for i = 1:n
    for j = i+1:n
        if B1(i, j) ~= 0
            nonZeroIndices = [nonZeroIndices; i, j];
        else
            zeroIndices = [zeroIndices; i, j];
        end
    end
end
% Shuffle the non-zero and zero indices randomly
shuffledNonZeroIndices = nonZeroIndices(randperm(size(nonZeroIndices, 1)), :);
shuffledZeroIndices = zeroIndices(randperm(size(zeroIndices, 1)), :);
% Choose the first two non-zero and zero indices from the shuffled arrays
numChosen = 2;
chosenNonZeroIndices = shuffledNonZeroIndices(1:numChosen, :);
chosenZeroIndices = shuffledZeroIndices(1:numChosen, :);
% % Define values to replace the chosen positions
% replacementValueNonZero = 999;  % Replace chosen non-zero positions with this value
% replacementValueZero = -999;     % Replace chosen zero positions with this value
% Replace the chosen non-zero positions with the replacement value
for i = 1:2
    B1(chosenNonZeroIndices(i, 1), chosenNonZeroIndices(i, 2)) = 0;
end
% Replace the chosen zero positions with the replacement value
for i = 1:2
    B1(chosenZeroIndices(i, 1), chosenZeroIndices(i, 2)) = 1;
end
% Display the modified matrix B
disp("Modified Matrix B1:");
disp(A);
M=B1+B1';
disp(M)
%G=graph(M)
%total=nnz(M);
%plot(G)
%%
% function

x0 = rand(20,1);
tspan=[0 1000];
[t,x] = ode45(@(t,x)fun(t,x,K),tspan, x0);

function xprime = fun(t, x, K)
    theta = 0.1;dh = 0.2;phi = 0.3;ita = 0.15;
    nn=20;
    
    i = 1; f = [];
    while i< nn  
        f = [f; (x(i)*(1 - theta*x(i))  - (x(i)*x(i+1))/(1 + x(i)));((phi*x(i)*x(i+1))/(1 + x(i)) - ita*x(i+1))];
        i = i + 1;  
    
    end
    xprime = f+dh*K*x;  
end
%%

.

[Rahim]

sorry here matrix M, not K ,I mean M=K.But sir, still the code is not running .Please help me

n = 20;

C= zeros(n, n);

C(1,n)=1;

C(1,n-1)=1;

C(2,n)=1;

A = zeros(n, n);

for i = 1:n

for j = 1:n

if j==i+1 || j == i + 2

A(i,j) = 1;

else

A(i, j) = 0;

end

end

end

B1=A+C;

B=B1+B1';

% Find non-zero and zero indices in the upper triangular part of the matrix

nonZeroIndices = []; % Initialize an array to store non-zero indices

zeroIndices = []; % Initialize an array to store zero indices

for i = 1:n

for j = i+1:n

if B1(i, j) ~= 0

nonZeroIndices = [nonZeroIndices; i, j];

else

zeroIndices = [zeroIndices; i, j];

end

end

end

% Shuffle the non-zero and zero indices randomly

shuffledNonZeroIndices = nonZeroIndices(randperm(size(nonZeroIndices, 1)), :);

shuffledZeroIndices = zeroIndices(randperm(size(zeroIndices, 1)), :);

% Choose the first two non-zero and zero indices from the shuffled arrays

numChosen = 2;

chosenNonZeroIndices = shuffledNonZeroIndices(1:numChosen, :);

chosenZeroIndices = shuffledZeroIndices(1:numChosen, :);

% % Define values to replace the chosen positions

% replacementValueNonZero = 999; % Replace chosen non-zero positions with this value

% replacementValueZero = -999; % Replace chosen zero positions with this value

% Replace the chosen non-zero positions with the replacement value

for i = 1:2

B1(chosenNonZeroIndices(i, 1), chosenNonZeroIndices(i, 2)) = 0;

end

% Replace the chosen zero positions with the replacement value

for i = 1:2

B1(chosenZeroIndices(i, 1), chosenZeroIndices(i, 2)) = 1;

end

% Display the modified matrix B

disp("Modified Matrix B1:");

disp(A);

M=B1+B1';

disp(M)

%G=graph(M)

%total=nnz(M);

%plot(G)

%%

% function

x0 = rand(20,1);

tspan=[0 1000];

[t,x] = ode45(@(t,x)fun(t,x,M),tspan, x0);

function xprime = fun(t, x, M)

theta = 0.1;dh = 0.2;phi = 0.3;ita = 0.15;

nn=20;

i = 1; f = [];

while i< nn

f = [f; (x(i)*(1 - theta*x(i)) - (x(i)*x(i+1))/(1 + x(i)));((phi*x(i)*x(i+1))/(1 + x(i)) - ita*x(i+1))];

i = i + 1;

end

xprime = f+dh*M*x;

end

Answer 1

[Sam Chak]

Hi @Rahim

The code runs now by fixing nn, but the system stability depends on the randomly chosen initial values x0.

Update: The code is generalized for the number n, which must satisfy two conditions. Firstly, it must be an even number, meaning it is divisible by 2 without leaving a remainder. Secondly, it must be greater than or equal to 6. In other words, the number must fulfill both the requirement of being an even number and the condition of having a value of 6 or higher.

n        = 20;  % must be both an even number and have a value of 6 or higher
C        = zeros(n, n);
C(1,n)   = 1;
C(1,n-1) = 1;
C(2,n)   = 1;
A        = zeros(n, n);

for i = 1:n
    for j = 1:n
        if j == i+1 || j == i + 2
            A(i,j) = 1;
        else
            A(i, j) = 0;
        end
    end
end

B1 = A  + C;
B  = B1 + B1';

% Find non-zero and zero indices in the upper triangular part of the matrix
nonZeroIndices = [];    % Initialize an array to store non-zero indices
zeroIndices    = [];	% Initialize an array to store zero indices

for i = 1:n
    for j = i+1:n
        if B1(i, j) ~= 0
            nonZeroIndices = [nonZeroIndices; i, j];
        else
            zeroIndices = [zeroIndices; i, j];
        end
    end
end

% Shuffle the non-zero and zero indices randomly
shuffledNonZeroIndices = nonZeroIndices(randperm(size(nonZeroIndices, 1)), :);
shuffledZeroIndices    = zeroIndices(randperm(size(zeroIndices, 1)), :);

% Choose the first two non-zero and zero indices from the shuffled arrays
numChosen            = 2;
chosenNonZeroIndices = shuffledNonZeroIndices(1:numChosen, :);
chosenZeroIndices    = shuffledZeroIndices(1:numChosen, :);

% % Define values to replace the chosen positions
% replacementValueNonZero = 999;  % Replace chosen non-zero positions with this value
% replacementValueZero = -999;     % Replace chosen zero positions with this value
% Replace the chosen non-zero positions with the replacement value

for i = 1:2
    B1(chosenNonZeroIndices(i, 1), chosenNonZeroIndices(i, 2)) = 0;
end

% Replace the chosen zero positions with the replacement value
for i = 1:2
    B1(chosenZeroIndices(i, 1), chosenZeroIndices(i, 2)) = 1;
end

% Display the modified matrix B
% disp("Modified Matrix B1:");
% disp(A);
M = B1 + B1';
% disp(M)

%%
tspan  = [0 4];
x0     = rand(n, 1);                                % <-- change in this line
[t, x] = ode45(@(t, x) fun(t, x, M, n), tspan, x0); % <-- change in this line
plot(t, x), grid on
xlabel('t'), ylabel('x(t)')

function xprime = fun(t, x, M, n)   % <-- change in this line
    theta  = 0.1;
    dh     = 0.2;
    phi    = 0.3;
    ita    = 0.15;
    nn     = 1/2*n + 1;             % <-- change in this line
    i      = 1; 
    f      = [];
    while i < nn  
        f  = [f;
              x(i)*(1 - theta*x(i))  - (x(i)*x(i+1))/(1 + x(i));
              (phi*x(i)*x(i+1))/(1 + x(i)) - ita*x(i+1)];
        i  = i + 1;  
    end
    xprime = f + dh*M*x;
end

[Rahim]

sir, x(t) and f both are 20*1 column vector and here n=20;and M is 20*20 adjacency matrix i want to solve this system .If you need more information please tell me I am here actively I will cooperate .I really appreciated for you help

Thank you so much sir ,

[Sam Chak]

Hi @Rahim, Then, just assign n = 20, as in the original code. The rest of the code is unaffected. Please test the code if it works in your machine.

[Rahim]

Sir, If I take time span[0 100],then solution should be x(t)=20*100 matrix write ?but here I am not getting this

please check again

[Sam Chak]

Hi @Rahim, the x(t) returned by ode45(), is an array, where each row in x corresponds to the solution at the value returned in the corresponding row of t. Click this icon and paste your code in the grey field, followed by clicking this icon to run the MATLAB code. We will check if there is any error.

[Rahim]

Questions1: Sir,if my function if the form .is it the right way to right? dp,dh are parameters.

Questions2: my second questtions is ,when we run this code sometimes x(t) is 145*20,sometimes x(t) is 165*20 why ?

function xprime = fun(t, x, M, n)

theta = 0.1;

dh = 0.2;

dp=0.25; %<......i have changed this

phi = 0.3;

ita = 0.15;

nn = 20; i = 1;

f = [];

while i < nn

f = [f;

x(i)*(1 - theta*x(i)) - (x(i)*x(i+1))/(1 + x(i))+dh*M*x; %<......i have changed this

(phi*x(i)*x(i+1))/(1 + x(i)) - ita*x(i+1)+dp*M*x]; %<......i have changed this

i = i + 1;

end

xprime = f ;

end

[Sam Chak]

Hi @Rahim

Reply to Q1. The setting "nn = 20" will not work if it is a 20th-order system. So, do not change the formula for nn in the code. Only set n = 20 at the beginning of the script. The choices of the parameters dh and dp are up to you. I have no mathematical and scientific knowledge about the system.

Reply to Q2. The ode45() algorithm uses the adaptive step size Runge-Kutta 4 method. Thus, randomized initial values result in a varying step size value in each run. If you want to fix the number of elements, refer to the "ODE solver" section.

n        = 20;                  % <-- fix at this value

⋮

%% ODE solver
tspan  = linspace(0, 4, 401);   % <-- set the number of elements from 0 to 4 sec
x0     = rand(n, 1);
[t, x] = ode45(@(t, x) fun(t, x, M, n), tspan, x0); 
plot(t, x), grid on
xlabel('t'), ylabel('x(t)')

% check number of elements in the state vector, x(t)
size(x)

ans = 1×2

401 20

%% The dynamic system
function xprime = fun(t, x, M, n)
    theta  = 0.1;
    dh     = 0.2;
    phi    = 0.3;
    ita    = 0.15;
    nn     = 1/2*n + 1;         % <-- Do not change this formula
    i      = 1; 
    f      = [];
    while i < nn  
        f  = [f;
              x(i)*(1 - theta*x(i))  - (x(i)*x(i+1))/(1 + x(i));
              (phi*x(i)*x(i+1))/(1 + x(i)) - ita*x(i+1)];
        i  = i + 1;  
    end
    xprime = f + dh*M*x;
end

[Rahim]

thanks