Getting a system of equations and the using ODE45

Hi, I have a coefficient matrix -

Xi=
   1.0e+03 *

    0.3002    0.3832    0.0577
         0         0         0
         0         0         0
         0         0         0
   -0.1425   -0.1576   -0.0162
         0         0         0
         0         0         0
   -0.0352   -0.0601   -0.0140
         0         0         0
         0         0         0
   -0.0941   -0.1304   -0.0245
    0.1394    0.1897    0.0296
    0.7197    0.8981    0.1662
    0.0058    0.0024    0.0061
   -1.0007   -1.2808   -0.2519
         0         0         0
   -0.0211   -0.0250   -0.0060
    0.3760    0.4933    0.0964
   -0.0242   -0.0253    0.0097
         0         0         0
% 

To get the system of ODEs I first defined the variables, and then multiplied it by transpose(Xi) to get the system-

syms x y z xx xy xz yy yz zz xxx xxy xxz xyy xyz xzz yyy yyz yzz zzz
var = [1;x;y;z;xx;xy;xz;yy;yz;zz;xxx;xxy;xxz;xyy;xyz;xzz;yyy;yyz;yzz;zzz]
A = transpose(Xi)
Xdot = A*var

where Xdot is the column containing the derivatives, ie. Xdot = [x' ; y' ; z']. The code above gives me the output:

Xdot =
 
      (4905784495475543*xxy)/35184372088832 - (6622198076817333*xxx)/70368744177664 - (1253707074885461*xx)/8796093022208 + (6330728071069091*xxz)/8796093022208 + (6526091537672905*xyy)/1125899906842624 - (550152702477171*xyz)/549755813888 - (4953639297776579*yy)/140737488355328 - (743004146196139*yyy)/35184372088832 + (6615295418414269*yyz)/17592186044416 - (3404894971584935*yzz)/140737488355328 + 2640157789935161/8796093022208
   (6675905246377611*xxy)/35184372088832 - (4586552884284477*xxx)/35184372088832 - (2773298922331439*xx)/17592186044416 + (987449040685507*xxz)/1099511627776 + (2666532593051851*xyy)/1125899906842624 - (5632811008410909*xyz)/4398046511104 - (8464188990339351*yy)/140737488355328 - (7048704279715849*yyy)/281474976710656 + (4339015521236507*yyz)/8796093022208 - (7129115399973291*yzz)/281474976710656 + 3370686532933275/8796093022208
(2082177599331991*xxy)/70368744177664 - (863749208859729*xxx)/35184372088832 - (4564685191010399*xx)/281474976710656 + (5847277141854391*xxz)/35184372088832 + (6877937169110265*xyy)/1125899906842624 - (8862390418139439*xyz)/35184372088832 - (491773253308733*yy)/35184372088832 - (3353395911004425*yyy)/562949953421312 + (6782944326528325*yyz)/70368744177664 + (5481194070111587*yzz)/562949953421312 + 2029864570192639/35184372088832
 

My question is

a) Is there a better way to define the system of ODEs?

b) How do I solve this system using ODE45? I do not want to manually input the equations, as my coefficient matrix will vary if I have a different dataset, so I cannot manually input equations everytime.

c) If I have 5 variables, say, x1,x2,x3,x4,x5 can I still use ODE45 to solve a system of ODEs containing 5 equations?

[Torsten]

As far as I can see, you have 3 equations for 19 unknowns. What are the missing 16 equations ?

[Shreshtha Chaturvedi]

It's a non-linear system of 3 equations and 3 variables/unknowns. xx=x^2 and so on.

Answer 1

[Torsten]

fun = @(x) Xi.'*[1;x(1);x(2);x(3);x(1)^2;x(1)*x(2);x(1)*x(3);x(2)^2;x(2)*x(3);x(3)^2;x(1)^3;x(1)^2*x(2);x(1)^2*x(3);x(1)*x(2)^2;x(1)*x(2)*x(3);x(1)*x(3)^2;x(2)^3;x(2)^2*x(3);x(2)*x(3)^2;x(3)^3];
tspan = [0 1];    % time interval of integration
x0 = [1; 1; 1];   % Initial conditions
[T,X] = ode45(fun,tspan,x0)  % solver call

[Shreshtha Chaturvedi]

Thank you, I tried it but it's giving me the following error-

Too many input arguments.

Error in odearguments (line 92)
f0 = ode(t0,y0,args{:});   % ODE15I sets args{1} to yp0.

Error in ode45 (line 107)
  odearguments(odeIsFuncHandle,odeTreatAsMFile, solver_name, ode, tspan, y0, options, varargin);

Error in TFR (line 35)
[T,X] = ode45(fun,tspan,x0)  % solver call
 
>> 

[Torsten]

Xi= 1.0e+03 * ...
    [0.3002    0.3832    0.0577
    0         0         0
    0         0         0
    0         0         0
    -0.1425   -0.1576   -0.0162
    0         0         0
    0         0         0
    -0.0352   -0.0601   -0.0140
    0         0         0
    0         0         0
    -0.0941   -0.1304   -0.0245
    0.1394    0.1897    0.0296
    0.7197    0.8981    0.1662
    0.0058    0.0024    0.0061
    -1.0007   -1.2808   -0.2519
    0         0         0
    -0.0211   -0.0250   -0.0060
    0.3760    0.4933    0.0964
    -0.0242   -0.0253    0.0097
    0         0         0];
fun = @(t,x) Xi.'*[1;x(1);x(2);x(3);x(1)^2;x(1)*x(2);x(1)*x(3);x(2)^2;x(2)*x(3);x(3)^2;x(1)^3;x(1)^2*x(2);x(1)^2*x(3);x(1)*x(2)^2;x(1)*x(2)*x(3);x(1)*x(3)^2;x(2)^3;x(2)^2*x(3);x(2)*x(3)^2;x(3)^3];
tspan = [0 1];    % time interval of integration
x0 = [1; 1; 1];   % Initial conditions
[T,X] = ode15s(fun,tspan,x0);  % solver call

Warning: Failure at t=1.340286e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.775558e-17) at time t.

plot(T,X(:,1))

[Walter Roberson]

Xi = 1.0e+03 * [
    0.3002    0.3832    0.0577
         0         0         0
         0         0         0
         0         0         0
   -0.1425   -0.1576   -0.0162
         0         0         0
         0         0         0
   -0.0352   -0.0601   -0.0140
         0         0         0
         0         0         0
   -0.0941   -0.1304   -0.0245
    0.1394    0.1897    0.0296
    0.7197    0.8981    0.1662
    0.0058    0.0024    0.0061
   -1.0007   -1.2808   -0.2519
         0         0         0
   -0.0211   -0.0250   -0.0060
    0.3760    0.4933    0.0964
   -0.0242   -0.0253    0.0097
         0         0         0];

% 
fun = @(t, x) Xi.'*[1;x(1);x(2);x(3);x(1)^2;x(1)*x(2);x(1)*x(3);x(2)^2;x(2)*x(3);x(3)^2;x(1)^3;x(1)^2*x(2);x(1)^2*x(3);x(1)*x(2)^2;x(1)*x(2)*x(3);x(1)*x(3)^2;x(2)^3;x(2)^2*x(3);x(2)*x(3)^2;x(3)^3];
tspan = [0 1];    % time interval of integration
x0 = [1; 1; 1];   % Initial conditions
[T,X] = ode45(fun,tspan,x0);  % solver call

Warning: Failure at t=1.349254e-02. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (2.775558e-17) at time t.

for K = 1 : size(X,2)
    figure(); plot(T, X(:,K)); title("X" + K);
end

[Shreshtha Chaturvedi]

I tried doing the same for 5 variables:

fun = @(t,x) Xi.'*[1;x(1);x(2);x(3);x(4);x(5);x(1)^2;x(1)*x(2);x(1)*x(3);x(1)*x(4);x(1)*x(5);x(2)^2;x(2)*x(3);x(2)*x(4);x(2)*x(5);x(3)^2;x(3)*x(4);x(3)*x(5);x(4)^2;x(4)*x(5);x(5)^2;x(1)^3;x(1)^2*x(2);x(1)^2*x(3);x(1)^2*x(4);x(1)^2*x(5);x(1)*x(2)^2;x(1)*x(2)*x(3);x(1)*x(2)*x(4);x(1)*x(2)*x(5);x(1)*x(3)^2;x(1)*x(3)*x(4);x(1)*x(3)*x(5);x(1)*x(4)^2;x(1)*x(4)*x(5);x(1)*x(5)^2;x(2)^3;x(2)^2*x(3);x(2)^3;x(2)^2*x(3);x(2)^2*x(4);x(2)^2*x(5);x(2)*x(3)^2;x(2)*x(3)*x(4);x(2)*x(3)*x(5);x(2)*x(4)^2;x(2)*x(4)*x(5);x(2)*x(5)^2;x(3)^3;x(3)^2*x(4);x(3)^2*x(5);x(3)*x(4)^2;x(3)*x(4)*x(5);x(3)*x(5)^2;x(4)^3;x(4)^2*x(5);x(4)*x(5)^2;x(5)^3];
tspan = [1 10000];    % time interval of integration
x0 = [249.270160982681; 80.9237201495384; 74.2320286869611; 79.7728847180895; 255.583672846571];   % Initial conditions
[T,X] = ode45(fun,tspan,x0)  % solver call

where Xi is a 56*5 matrix. But when I try to run the last line (ode45) it says

Error using  * 
Incorrect dimensions for matrix multiplication. Check that the number of columns in the first
matrix matches the number of rows in the second matrix. To operate on each element of the matrix
individually, use TIMES (.*) for elementwise multiplication.

However, the problem is not with the matrix multiplication as if I comment the last line and run it, my program executes without any error. Any hints as to why is it behaving this way when I change dimensions?

[Torsten]

Your vector is 58x1 instead of 56x1.

Answer 2

[Walter Roberson]

I recommend that you read the first example for odeFunction() to see the flow to turn an array of equations into something that can be evaluated numerically by ode45()