Numerical Methods for Singular systems

Could you please take a look at the code I wrote for the delayed singular system class and let me know if it's correct? Thank you all in advance!

%% Phase plane
clear all
clc

%% conditions


A_x=[0.5 0.7 -0.5;0.8 1.5 0.2; 0.9 0.8 0.3];
A_d=[1.6 0.58 0.4;0.96 1.7 0.8;0.5 0.8 0.45];
B=[1.4;0.5; 0.7];


stime=0;
endtime=600;
h=0.01;     
t=-stime:h:endtime;
N_0=stime/h;
N_1=endtime/h;
N=N_0+N_1;
%bound=110;




for i=1:N+1 
          
    
          
    if i<=N_0+1;
           
        x1(:,i)=0; 
        x2(:,i)=0;  
        x3(:,i)=0;   
    else
     d_1=floor(1.5+sin(0.5*i*h));

    
  if i<=3
      
     x1(:,i)=x1(:,i-1)+h*(A_x(1,1)*x1(:,i-1)+A_x(1,2)*x2(:,i-1)+A_x(1,3)*x3(:,i-1)+B(1,1)*sin(0.5*i));
               
     x2(:,i)=x2(:,i-1)+h*(A_x(2,1)*x1(:,i-1)+A_x(2,2)*x2(:,i-1)+A_x(2,3)*x3(:,i-1)+B(2,1)*sin(0.5*i));

     x3(:,i)=-inv(A_x(3,3))*(A_x(3,1)*x1(:,i)+A_x(3,2)*x2(:,i)+B(3,1)*sin(0.5*i));
          
  else
       x1(:,i)=x1(:,i-1)+h*(A_x(1,1)*x1(:,i-1)+A_x(1,2)*x2(:,i-1)+A_x(1,3)*x3(:,i-1)+A_d(1,1)*x1(:,i-d_1-1)+A_d(1,2)*x2(:,i-d_1-1)+A_d(1,3)*x3(:,i-d_1-1)+B(1,1)*sin(0.5*i));
               
       x2(:,i)=x2(:,i-1)+h*(A_x(2,1)*x1(:,i-1)+A_x(2,2)*x2(:,i-1)+A_x(2,3)*x3(:,i-1)+A_d(2,1)*x1(:,i-d_1-1)+A_d(2,2)*x2(:,i-d_1-1)+A_d(2,3)*x3(:,i-d_1-1)+B(2,1)*sin(0.5*i));
       
       x3(:,i)=-inv(A_x(3,3))*(A_x(3,1)*x1(:,i)+A_x(3,2)*x2(:,i)+A_d(3,1)*x1(:,i-d_1)+A_d(3,2)*x2(:,i-d_1)+A_d(3,3)*x3(:,i-d_1)+B(3,1)*sin(0.5*i));

  end           
    end
  
end

%% Plotting Graphs

figure(1)
% grid on 
% hold on
% box on
plot(t,x1(1,:),'b-','linewidth',1)
hold on
plot(t,x2(1,:),'r:','linewidth',2)
plot(t,x3(1,:),'r:','linewidth',2)
 xlabel('Time (t)')
%ylabel('Responses x(t)')
 legend('x_{1}(t)','x_{2}(t)');
% xlim([0,1000]);
% %ylim([-0.2,0.8]);

%

Answer 1

[Torsten]

The loop index i does not equal time t. So expressions like

d_1=floor(1.5+sin(0.5*i*h));
x1(:,i-d_1-1)
x2(:,i-d_1-1)
x3(:,i-d_1-1)
sin(0.5*i)

in your code are clearly wrong.

Up to t ~ 2.43, the delay part is zero, and the results should be the same as with the code below. It can serve as validation code for your results after you've made the necessary corrections.

M = [1 0 0;0 1 0;0 0 0];
A = [0.5 0.7 -0.5;0.8 1.5 0.2; 0.9 0.8 0.3];
B = [1.4;0.5; 0.7];
tstart = 0;
tend = 2.43;
tspan = [tstart,tend];
f = @(t,y)A*y+B*sin(0.5*t);
y0 = [0;0;0];
options = odeset('Mass',M);
[T,Y] = ode15s(f,tspan,y0,options);
plot(T,Y)
grid on

●

[Sam Chak]

Hi @Torsten and @Walter Roberson

I did not use the original code provided by the OP. Instead, I am writing from scratch. But I am unsure if the following approach is the correct method for solving delay differential algebraic equations in MATLAB. Could you please review it?

For differential-algebraic equations (DAEs), we typically use the ode15s or ode23t solvers, and for delay differential equations (DDEs), we commonly utilize either the dde23 or ddesd solvers. However, for a combination of both types, I am uncertain which solver to choose.

tspan   = [0 5];
M       = [1 0 0    % singular mass matrix
           0 1 0
           0 0 0];
options = odeset('Mass', M);
sol     = ddesd(@ddefun, @delx, @history, tspan, options);

%% plot results
plot(sol.x, sol.y, '-o'), grid on
legend('x_1', 'x_2', 'x_3', 'location', 'northwest')
title('Solution for delay differential algebraic system')
xlabel('Time t')
ylabel('Solution x')

%% delay function
function d = delx(t, x)
    d = t - 1.5 - sin(0.5*t);
end

%% delay differential equations
function dydt = ddefun(t,x,Z)
    A    = [0.5,  0.7, -0.5
            0.8,  1.5,  0.2
            0.9,  0.8,  0.3];
    
    Ad   = [1.60, 0.58, 0.40
            0.96, 1.70, 0.80
            0.50, 0.80, 0.45];
    
    B    = [1.4
            0.5
            0.7];
    
    dydt = A*x + Ad*Z + B*sin(0.5*t);
end

%% solution history 
function s = history(t)
    s = zeros(3, 1);
end

[Torsten]

@Sam Chak

Your code would be correct of "ddesd" worked for differential-algebraic delay differential equations. But it only works for differential delay equations as far as I know. Thus the "0" in front of x3dot is the problem.

ddeset

          AbsTol: [ positive scalar or vector {1e-6} ]
          Events: [ function_handle ]
     InitialStep: [ positive scalar ]
        InitialY: [ vector ]
           Jumps: [ vector ]
         MaxStep: [ positive scalar ]
     NormControl: [ on | {off} ]
       OutputFcn: [ function_handle ]
       OutputSel: [ vector of integers ]
          Refine: [ positive integer {1} ]
          RelTol: [ positive scalar {1e-3} ]
           Stats: [ on | {off} ]

[Le]

Ok, Thank you very much @Torsten

[Torsten]

Could you find the solution for the system without the delay terms (thus up to 2.4389) (e.g. using the implicit Euler method) ?

f = @(x) x-1.5-sin(0.5*x);
fzero(f,2)

ans = 2.4389

[Le]

Could you please help me correct the code above @Torsten? Thank you very much."

[Torsten]

You post your questions about code for delay differential equations for almost one year now, but it seems you don't make any progress (at least in the programming part). It's not possible to give solutions to obvious homework problems, but you should enhance your programming skills: try MATLAB Onramp - an introduction to MATLAB free of costs to learn the basics of the new language:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted

[Le]

ok, thank you @Torsten so much!

Answer 2

[Walter Roberson]

    if i<=N_0+1;
           
        x1(:,i)=0; 
        x2(:,i)=0;  
        x3(:,i)=0;   
    else
     d_1=floor(1.5+sin(0.5*i*h));
    
  if i<=3

This smells like you forgot the "end" for the first "if". In MATLAB, nesting is not determined by indentation.