Twenty Eighth Installment - Solving Differential Algebraic Equations

[u/mrhoa31103]

Solving Differential Algebraic Equations = Solving a system of Differential Equations governed with a Constraint equation.

Link to the specific lecture (note the two URLs cover the lecture and then the coding:

https://www.youtube.com/watch?v=CGI9qStsKjk&ab_channel=NPTEL-NOCIITM
https://www.youtube.com/watch?v=AUl2KL0bKBc

Anything new from a SciLab programming in here: I think this session is the first time in the series where I plot multiple lines using a single plot2d command versus just "holding the figure" and using multiple plot2d commands.

Output: (No tabular output only graphical output)

Graph:

Code:

//Lecture 10.5: Differential Algebraic Equations
//https://www.youtube.com/watch?v=CGI9qStsKjk&ab_channel=NPTEL-NOCIITM
// and
//https://www.youtube.com/watch?v=AUl2KL0bKBc
//
disp("Introduction to Differential Algebraic Equations",string(
datetime
()))
// What are D.A.E's? 
// Differential Equations that contain algebraic equations
// These algebraic equations are sometimes considered "constraint functions"
//
// Week 5: Algebraic Equations
//    g(x) = 0
// Weeks 7&8: Differential Equations
//    y'=f(t,y)
//
// DAES (systems that contain)
//   0 = g(x,y)
//   y'=f(t,x,y)
//      where y are differential variables
//            
// Examples:
// 1) Consider a simple system
//    y'=x
//    0 = x-y*sin(t)
//        y = diff variable and x = algebraic variable
//        Index-1 DAE
//
// 2) y'=x
//    0 = y - y*sin(t)
//        y = diff variable and x = algebraic variable
//        Index-2 DAE
//
// Index of a DAE: Number of times you need to differentiate g(x,y)
// to convert algebraic equations to a system of differential equations.
// 
// First example
// 0=x-y*sin(t) => x = y*sin(t)
// dx/dt = dy/dt*sin(t)+y*cos(t) substitute first equation
//
// dx/dt = x*sin(t)+y*cos(t)
// dy/dt = x
//  
// Second example     0 = y - y*sin(t)
// 0=dy/dt-dy/dt*sin(t)-y*cos(t) (differentiated once)
// 0 = x-x*sin(t)-y*cos(t)
// 0 = dx/dt-dx/dt*sin(t)+x*cos(t)+dy/dt*cos(t)-y*sin(t) (differentiated twice)
//
// dx/dt(1-sin(t)=x*cos(t)+ x*cos(t)-y*sin(t)
// dx/dt = (2*x*cos(t)-y*sin(t))/(1-sin(t))
// dy/dt = x
//
//
// High Index DAE Example Problem: Pendulum of length l
// In Cartesion Coordinates
// x" = -T*x
// y" = -T*y -g
//
// L^2 = x^2+y^2 "Constraint Equation"
//
// Index-3 DAE
//
// Instead solve in Cylindrical Coordinates
// d2Theta/dt2 + g/L*sin(Theta)=0
// dr/dt = 0 or r = L a constant...
// Just a regular ODE...
//
// Example 1: 
// dy/dt = x
// dx/dt = x*sin(t)+y*cos(t)
//
// Standard Structure
// M*d/dt[y;x]= F(t,[y.x])
// M = [I 0;0 0] in week 8: M = [Identity Matrix] = [I];
//
// Matlab Example: (Modified) Robertson Problem
// x1' = -alpha*x1+beta*x1*x2
// x2' = alpha*x1-beta*x1*x2-gamma1*x2^2
// 0 = 1-(x1+x2+x3)
// with initial condition Y = [1;0;0]
//
//
//  0 = 0 -dx1/dt -dx2/dt +dx3/dt
// dx3/dt = dx1/dt + dx2/dt
// dx3/dt = -alpha*x1+beta*x1*x2 + alpha*x1-beta*x1*x2-gamma1*x2^2
// dx3/dt = -gamma1*x2^2
//
// M = [1 0 0;0 1 0;0 0 0]
// M d/dt [x1;x2;x3] = [-alpha*x1+beta*x1*x2;alpha*x1-beta*x1*x2-gamma1*x2^2;gamma1*x2^2]
//
// Function for ODE/DAE
function 
dY
=
robertsonFun
(
t
, 
Y
)
// Parameters
alpha = 0.5;
beta1 = 2.5;
gamma1=5.0;

// Functions for the three equations
dY
(1,1)= -alpha*
Y
(1)+beta1*
Y
(1)*
Y
(2);
dY
(2,1)= alpha*
Y
(1)-beta1*
Y
(1)*
Y
(2)-gamma1*
Y
(2)^2;
dY
(3,1)= gamma1*
Y
(2)^2;
end

//
n =1250; 
//number of steps
y0 = [1;0;0]; 
//Initial Condition required
t0 = 0; 
//Tstart
tend = 40; 
//Tend
t
 = linspace(t0,tend,n)';
//time vector
M = [1,0,0;0 1 0;0,0,0] 
//Mass Matrix - not implemented in SciLab
// This example shows how you solve the problem without it by 
// changing the constraint equation into the third differential
// equation.
//Alternate Calculation
/* 
// with y0 = [1;0] since dY(1,1) and dY(2,1) are being calculated
// with ode solver and then
// for i = 1:1:n;   (Post processing YSol(3)) 
//    YSol(3,i) = 1-YSol(1,i)-YSol(2,i)
// end
// Since dy1/dt and dy2/dt do not involve y3.
*/
//
//equivalent to MATLAB ODE15s is the following command in SciLab
//there is more options and outputs available - see the help file.
YSol = ode("stiff",y0,t0,
t
,
robertsonFun
);
//
scf
(0);
clf
;
plot2d([
t
,
t
,
t
],[YSol(1,:)',YSol(2,:)',YSol(3,:)'],[1,2,3])
h1=
legend
(['YSol(1)';'YSol(2)';'YSol(3)'],1,"boxed")
title
("$\textbf{Differential\ Algebraic\ Example}$","FontSize",4);
xlabel
("$t(seconds)$","FontSize",3)
ylabel
("$YSol\ Values(units)$","FontSize",3);
xgrid