ode45 : Unable to perform assignment because brace indexing is not supported for variables of this type

I am still not able to understand what I am doing wrong !!

My problem is to make 2 plots x1(t) and x2(t) vs t, and x1(t) vs x2(t).

I have x1(0) = x2(0) = 1, my time span is 0 to 140 seconds and

My 2 non linear differential equations are:

My input u is:

My code is:

clc
clear all
close all

a = 2;
b = 3;
k = 1;
% input function
u = @(x1, x2) -(x1*tanh(x1) + k*x2 + k*a*x1 + k*x1^2 + (a^2)*(x1)*tanh(x1) + a*x2*tanh(x1) + 2*a*(x1^2)*tanh(x1) + 2*x1*x2*tanh(x1) + b*x1*x2)*(1+x2^2);
tspan = [0 140];
x0 = [1, 1]; % Initial Condition
[t,x] = ode45(@(x1, x2) odefcn(x1, x2, u, a, b), tspan, x0);

%Plots
figure
plot(t,x(:,1),'-r');
grid on;
title('State Response x1(t) and x2(t) w.r.t time');
xaxis('time(s)');
yaxis('x1(t) and x2(t) in m');
hold on;
plot(t,x(:,2),'-b');
legend('x1(t)', 'x2(t)')
hold off

figure
plot(x(:,1), x(:,2), '-b');
grid on;
title('State Response x1(t) vs x2(t)');
xaxis('x1(t) in m');
yaxis('x2(t) in m');


% Function Handle
function dxdt = odefcn(x1, x2, u, a, b)
dxdt = zeros(2,1);
dxdt(1) = (a*x1 + x2)*tanh(x1);
dxdt(2) = b*x1*x2 + (1/(1 +x2^2))*u;
end

Error I am getting is:

Unable to perform assignment because the left and right sides have a different number of elements.

Error in NLST_HW5_3>odefcn (line 37)
dxdt(1) = (a*x1 + x2)*tanh(x1);

Error in NLST_HW5_3>@(x1,x2)odefcn(x1,x2,u,a,b)

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

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in NLST_HW5_3 (line 12)
[t,x] = ode45(@(x1, x2) odefcn(x1, x2, u, a, b), tspan, x0);

I do not understand why the left and right hand sides have different elements since I am just inputing single value into dxdt(1) and dxdt(2), there are no matrices or anything being used.

I would be grateful if anyone can explain me the reason for error and also help solve the problem !!

Answer 1

[Walter Roberson]

Your x0 is a vector of length 2. You have two boundary conditions. ode45 is going to pass a vector of that length, 2, as the second parameter to the function, x2. The first parameter will be scalar time, x1.

Your calculation involves a*x1 (remember x1 is time) + x2 (a vector) so it gives a vector result.

You could

ode45(@(t,x) odefcn(x(1), x(2), u, a, b), etc

However notice that your u is a function handle, and that you are multiplying by u. You cannot multiply by a function handle. You can multiply by u(x1, x2) though.

[Walter Roberson]

I was wrong, you should not be using 4 initial conditions for this. For that a, b, k then

    @(t,x)[tanh(x(1)).*(x(1).*2.0+x(2));
           x(1).*-2.0-x(2)-x(1).^2.*tanh(x(1)).*4.0-ix(1).*tanh(x(1)).*5.0-x(2).*tanh(x(1)).*2.0-x(1).^2-x(1).*x(2).*tanh(x(1)).*2.0]

Using the symbolic toolbox:

syms x1(t) x2(t) u(t)

a = 2;
b = 3;
k = 1;

u(t) = -(x1*tanh(x1) + k*x2 + k * a * x1 + k * x1^2 + a^2 * x1 * tanh(x1) + a * x2 * tanh(x1) + 2 * a * x1^2 * tanh(x1) + 2 * x1 * x2 * tanh(x1) + b * x1 * x2) * (1 + x2^2);
dx1 = diff(x1);
dx2 = diff(x2);
eqn1 = dx1 == (a * x1 + x2) * tanh(x1);
eqn2 = dx2 == b * x1 * x2 + 1/(1+x2^2) * u;
ic(1) = 1;  %x1
ic(2) = 1;  %x2

[eqs, vars] = reduceDifferentialOrder([eqn1, eqn2], [x1, x2]);
[M,F] = massMatrixForm(eqs,vars);
f = M\F;
odefun = odeFunction(f,vars);
disp(odefun)

ode45(odefun, [0 10], ic)

[Japnit Sethi]

This does plot the x1 and x2 vs t correctly !! Can you explain your code a little bit more since the set up is very different than what my original code is !! Specifically these lines !!

[eqs, vars] = reduceDifferentialOrder([eqn1, eqn2], [x1, x2]);
[M,F] = massMatrixForm(eqs,vars);
f = M\F;
odefun = odeFunction(f,vars);

Also how can i plot x1 vs x2 ?

I am also curious of what was I doing wrong with my own code and how could i have corrected it !!

[Walter Roberson]

reduceDifferentialOrder converts derivatives into separate auxillary variables, so that each resulting equation is individually first-order (but possibly involving variables that are standing in for higher derivatives.) reduceDifferentialOrder is specific to the Symbolic Toolbox.

The massMatrixForm and \ is a way of rewriting the equations in terms of the variables that are standing in for the derivatives. massMatrixForm is specific to the Symbolic Toolbox.

odeFunction coverts the resulting expressions involving functions into an anonymous numeric function involving input variables. It is similar to matlabFunction that is used to convert symbolic expressions into anonymous numeric functions, but has the extra step of converting the function references from x1(t) into a input parameter reference. odeFunction is specific to the Symbolic Toolbox.

[Japnit Sethi]

Thanks for the help !!

My final code was:

% Using Symbolic Toolbox
syms x1(t) x2(t) u(t)
a = 2;
b = 3;
k = 1;

% Input Function
u(t) = -(x1*tanh(x1) + k*x2 + k * a * x1 + k * x1^2 + a^2 * x1 * tanh(x1)...
    + a * x2 * tanh(x1) + 2 * a * x1^2 * tanh(x1) + 2 * x1 * x2 * tanh(x1) + b * x1 * x2) * (1 + x2^2);

% Setting up first order differential equations x1' and x2'
dx1 = diff(x1);
dx2 = diff(x2);
eqn1 = dx1 == (a * x1 + x2) * tanh(x1);
eqn2 = dx2 == b * x1 * x2 + 1/(1+x2^2) * u;

%Initial Conditions for x1 and x2
ic(1) = 1;  %x1
ic(2) = 1;  %x2

% Reducing a higher order differential equation to an equivalent first
% order differential equations
[eqs, vars] = reduceDifferentialOrder([eqn1, eqn2], [x1, x2]);

% Converting differential equation to mass matrix form
% Way of rewriting the equations in form of variables and the constants
% being inserted, Note 'f' gives the final x1' and x2' eq after computation
[M,F] = massMatrixForm(eqs,vars);
f = M\F;

% odeFunction is used to convert symbolic expressions in our case 'f' to
% function handle in our case 'odefun' to be used for ODE solvers
odefun = odeFunction(f,vars);
disp(odefun)

% Using ode45 in time period of 0 to 140 second and with initial conditions
% ic
[t,x] =  ode45(odefun, [0 140], ic);

%Plots
figure
plot(t,x(:,1),'-r');
grid on;
title('State Response x1(t) and x2(t) w.r.t time');
xlabel('time(s)');
ylabel('x1(t) and x2(t) in m');
hold on;
plot(t,x(:,2),'-b');
legend('x1(t)', 'x2(t)')
hold off

figure
plot(x(:,1), x(:,2), '-k');
grid on;
title('State Response x1(t) vs x2(t)');
xlabel('x1(t) in m');
ylabel('x2(t) in m');

I still do not understand the reasoning behind the Mass matrix form and wondering how I could have done this problem without using the symbolic toolbox !!

Answer 2

[Steven Lord]

Let's step through your code section by section.

clc
clear all
close all

Many people start their code with these three lines, but I recommend only using each one when you actually need to clear the command window, clear all variables and breakpoints, or to close all figure windows. The biggest "hammer" of all three of those is clear all, and one way to avoid needing to use it is to write functions instead of scripts. Each function has its own workspace that is set up when you call the function and destroyed when that function call finishes. Think of a function workspace as a blank piece of paper -- each time you start writing on a new piece of paper, you usually don't need to worry about what you wrote on the last piece of paper.

a = 2;
b = 3;
k = 1;
% input function
u = @(x1, x2) -(x1*tanh(x1) + k*x2 + k*a*x1 + k*x1^2 + (a^2)*(x1)*tanh(x1) + ...
    a*x2*tanh(x1) + 2*a*(x1^2)*tanh(x1) + 2*x1*x2*tanh(x1) + b*x1*x2)*(1+x2^2);

I just made this one span two lines so you could see the whole thing without scrolling in the Answers interface.

tspan = [0 140];
x0 = [1, 1]; % Initial Condition
[t,x] = ode45(@(x1, x2) odefcn(x1, x2, u, a, b), tspan, x0);

So odefcn gets called with five inputs.

The first is the time input (denoted in the ode45 documentation as t.) It is a scalar.

The second is the vector of states for your ODE. Since you have two equations (and two initial conditions) this vector will have two inputs

The third is the function handle u, which itself takes two inputs.

The fourth and fifth are scalars.

%Plots
figure
plot(t,x(:,1),'-r');
grid on;
title('State Response x1(t) and x2(t) w.r.t time');
xaxis('time(s)');
yaxis('x1(t) and x2(t) in m');

There are no functions xaxis and yaxis in MATLAB. You want the xlabel and ylabel functions instead.

hold on;
plot(t,x(:,2),'-b');
legend('x1(t)', 'x2(t)')
hold off
figure
plot(x(:,1), x(:,2), '-b');
grid on;
title('State Response x1(t) vs x2(t)');
xaxis('x1(t) in m');
yaxis('x2(t) in m');

% Function Handle
function dxdt = odefcn(x1, x2, u, a, b)
dxdt = zeros(2,1);
dxdt(1) = (a*x1 + x2)*tanh(x1);
dxdt(2) = b*x1*x2 + (1/(1 +x2^2))*u;
end

Remember what ode45 will pass into your odefcn as inputs. From above:

The first is the time input (denoted in the ode45 documentation as t.) It is a scalar.

The second is the vector of states for your ODE. Since you have two equations (and two initial conditions) this vector will have two inputs

The third is the function handle u, which itself takes two inputs.

The fourth and fifth are scalars.

So (a*x1 + x2)*tanh(x1) returns a (scalar * scalar + two-element-vector) + scalar or a two element vector. You can't fit a two element vector on the right side of the equation into one element of the vector dxdt.

You have similar problems when you try to define dxdt(2). b*x1*x2 is a two element vector and you want to use array division and exponentiation (./ and .^) instead of matrix division and exponentiation (/ and ^). See this documentation page for more information.

The line where you define dxdt(2) also has a problem in your use of u. You can't multiply a number (scalar, vector, matrix, or array) by a function handle. You can multiply a number by the value you get when you evaluate a function handle, however. So you probably want to multiply by u(x2(1), x2(2)) or something similar.

[Japnit Sethi]

Thanks for a detailed response !! I am still trying to understand the ode45 !!

I want to go step by step

I changed my code to:

a = 2;
b = 3;
k = 1;
% input function
% u = @(x1, x2) -(x1*tanh(x1) + k*x2 + k*a*x1 + k*x1^2 + (a^2)*(x1)*tanh(x1) + ...
%     a*x2*tanh(x1) + 2*a*(x1^2)*tanh(x1) + 2*x1*x2*tanh(x1) + b*x1*x2)*(1+x2^2);
u = @(x) -(x(1)*tanh(x(1)) + k*x(2) + k*a*x(1) + k*x(1)^2 + (a^2)*(x(1))*tanh(x(1)) + ...
    a*x(2)*tanh(x(1)) + 2*a*(x(1)^2)*tanh(x(1)) + 2*x(1)*x(2)*tanh(x(1)) + b*x(1)*x(2))*(1+x(2)^2);
tspan = [0 140];
x0 = [1, 1]; % Initial Condition
[t,x] = ode45(@(t, x) odefcn(t, x, u, a, b), tspan, x0);

I first of all changed x1 and x2 (which are 2 variables) to x(1) and x(2) and put it in a single 2x1 vector.

My question in this part is that "Do I need to change from matrix multiplication (*) to array multiplication (.*) in my function handle u ?"

%Plots
figure
plot(t,x(:,1),'-r');
grid on;
title('State Response x1(t) and x2(t) w.r.t time');
xlabel('time(s)');
ylabel('x1(t) and x2(t) in m');
hold on;
plot(t,x(:,2),'-b');
legend('x1(t)', 'x2(t)')
hold off

figure
plot(x(:,1), x(:,2), '-b');
grid on;
title('State Response x1(t) vs x2(t)');
xlabel('x1(t) in m');
ylabel('x2(t) in m');

xlabel and ylabel was a pure rookie mistake so thanks for pointing that out !! no question here right now!!

% Function Handle
function dxdt = odefcn(t, x, u, a, b)
dxdt = zeros(4,1);
dxdt(1:2,:) = (a*x(1) + x(2))*tanh(x(1));
dxdt(3:4,:) = b*x(1).*x(2) + (1./(1 + (x(2).^2))).*u(x);
end

I changed from matrix multiplication and division (* and /) to array multiplication and division (.* and ./) here and used u as a function handle and thus defined it as u(x) as suggested by you !!

Also my question from your explanation "So (a*x1 + x2)*tanh(x1) returns a (scalar * scalar + two-element-vector) + scalar or a two element vector. You can't fit a two element vector on the right side of the equation into one element of the vector dxdt."

Did you mean: " (Scalar * two-element-vector + two-element-vector)* scalar or a two element vector " ?

So when I run it I am getting an error !!

Error using odearguments (line 95)
@(T,X)ODEFCN(T,X,U,A,B) returns a vector of length 4, but the length of initial conditions vector is 2. The vector returned by
@(T,X)ODEFCN(T,X,U,A,B) and the initial conditions vector must have the same number of elements.

Error in ode45 (line 115)
  odearguments(FcnHandlesUsed, solver_name, ode, tspan, y0, options, varargin);

Error in NLST_HW5_3 (line 15)
[t,x] = ode45(@(t, x) odefcn(t, x, u, a, b), tspan, x0);

My question is that the length of initial condition is obviosly gonna be a vector of length 2 sicne the initial condition is x1 and x2 at time 0 is 1 each i.e. x1(0) = x2(0) = 1 !!

And how can i fix this !! Thanks

[Walter Roberson]

You have two differential equations. In nearly all cases, for each differential equation, you need a boundary condition for every order at least as high as the order of your equation. In your case your equation is in the first derivative and you need a boundary condition for the derivative and another for the non-derived function. If you happened to be working with the second derivative in an equation that did not mention the first derivative, you would still need a boundary condition for the second derivative, and the first derivative, and the un-derived function.

Thus you should have four boundary conditions: x1(t0), x1'(t0), x2(t0), x2'(t0) though not necessarily in that order.

[Japnit Sethi]

Oh Ok !! So I found the initial conditions of the first derivatives using initial conditions of x(1) and x(2) which are x1'(0) = 2.2848 and x2'(0) = -13.9007 and thus changed my x0 to

% x0 = [1,  2.2848, 1, -13.9007240274249];
x0 = [1,  1, 2.2848, -13.9007240274249];

I have tried both of the orders, eventhough I am now getting plots but there is still something fishy going on !!

My plots are :

[Walter Roberson]

Approximately this is more likely:

D(1,1) = y(2)
D(3,1) = y(4)
D(2,1) = equation that starts with a
D(4,1) = equation that starts with b

[Japnit Sethi]

What do you mean by the first 2 equations i.e. y(4) and y(2)

[Japnit Sethi]

If you meant

dxdt(1,:) = x(1);
dxdt(3,:) = x(2);
dxdt(2,:) = (a*x(1) + x(2)).*tanh(x(1));
dxdt(4,:) = b*x(1).*x(2) + (1./(1 + (x(2).^2))).*u(x);

That still outputs wrong plots

[Walter Roberson]

You are using the wrong x entries. Remember you no longer have only two of them. Boundary for x1(t) is in x(1), and for x2(t) is in x(3)

[Japnit Sethi]

Can you point it out in my code ? Thanks

This is my current code:

clc
clear all
close all

a = 2;
b = 3;
k = 1;
% input function
u = @(x) -(x(1).*tanh(x(1)) + k*x(2) + k*a*x(1) + k*x(1).^2 + (a^2)*(x(1)).*tanh(x(1)) + ...
    a*x(2).*tanh(x(1)) + 2*a*(x(1)^2).*tanh(x(1)) + 2*x(1).*x(2).*tanh(x(1)) + b*x(1).*x(2)).*(1+x(2)^2);
tspan = [0 140];

% Initial Condition
x0 = [1,  2.2848, 1, -13.9007240274249];
% x0 = [1,  1, 2.2848, -13.9007240274249];
[t,x] = ode45(@(t, x) odefcn(t, x, u, a, b), tspan, x0);

%Plots
figure
plot(t,x(:,1),'-r');
grid on;
title('State Response x1(t) and x2(t) w.r.t time');
xlabel('time(s)');
ylabel('x1(t) and x2(t) in m');
hold on;
plot(t,x(:,3),'-b');
legend('x1(t)', 'x2(t)')
hold off

figure
plot(x(:,1), x(:,3), '-k');
grid on;
title('State Response x1(t) vs x2(t)');
xlabel('x1(t) in m');
ylabel('x2(t) in m');


% Function Handle
function dxdt = odefcn(t, x, u, a, b)
dxdt = zeros(4,1);
dxdt(1,:) = x(1);
dxdt(3,:) = x(2);
dxdt(2,:) = (a*x(1) + x(2)).*tanh(x(1));
dxdt(4,:) = b*x(1).*x(2) + (1./(1 + (x(2).^2))).*u(x);
end

[Walter Roberson]

dxdt(4,:) = b*x(1).*x(2) + (1./(1 + (x(2).^2))).*u(x);

That line assumes that x(1) is the boundary condition for x1(t) and that x(2) is the boundary condition for x2(t) . But x(2) pertains to x1'(t) not to x2(t)