Solving ODEs using trapezoidal method in vector/matrix notation

1 view (last 30 days)
Brenton Foreman
Brenton Foreman on 22 Oct 2021
Answered: Chris on 22 Oct 2021
I have been given an example of solving multiple differential ODEs using the numerical trapezoidal method, and encountering issues when attempting to index vector and matrix values inside a for loop.
Here is the example:
My code is as follows:
clear; clc;
Tend = 0.1; %End time
N = 1000; %Number of steps
T = Tend/N; %step size
t = 0:T:Tend; %time span
A = [0 -5.38e6; 1 -47.106]; %System Matrix
B = [5.38e6; 0]; %Input Matrix
I = eye(2);
C1 = inv(I - A.*(T/2));
C2 = (I + A.*(T/2));
C3 = C1.*C2; %x(t) coefficient
C4 = (B.*(T/2)); %u(t) coefficient
% Initial Conditions
x1(1) = 0;
x2(1) = 0;
%Create Cell Array
x = zeros(2,1);
xx = size(x);
y = cell(xx);
%Functions
xdot1 = @(x1,t) x1;
xdot2 = @(x2,t) x2;
xdot = @(x1,x2,t)[xdot1;xdot2];
%Assign Values to Cell Array
y{1} = xdot1;
y{2} = xdot2;
%Input
ut = @(t) 187.8.*cos(377.*t);
%Pre-Allocate and store results
xt1 = zeros(2,N+1);
%Run and store values
for n=1:N
xt1(n+1)= C3.*[y{1}; y{2}] + C4.*(ut(t(n)+ut(t(n+1))));
end
I am recieving the following error:
Error using vertcat
Nonscalar arrays of function handles are not allowed; use cell arrays instead.
Any help at finding my mistake would be greatly appreciated

Answers (1)

Chris
Chris on 22 Oct 2021
A function handle works like this:
f = @(x,t) x.^2;
y = f(3,12);
the function f takes 3 and 12 as input variables and names them x and t, respectively. Since t is unused in the function, it's irrelevant and may not need to be in the function. It squares x and returns 9. equivalently, f = @(x) x.^2;
A couple things I see:
xdot1 = @(x1,t) x1;
You don't have a function x1 defined--x1 is a vector with an initial value. Calling xdot1(someX,someTime) would always return the entire vector x1 (initially 0). Similarly for xdot2.
xt1(n+1)= C3.*[y{1}; y{2}] + ...
y{1}, at this point, is a function
y{1} = @(x1,t) x1
It can't be concatenated into an array. If you give it numerical input, it will return a number, which is possible to use in an array.
xt1(n+1)= C3.*[y{1}(3,5); y{2})4,2)] + ...
Does this help? I'm sure there will be some more work to do once these errors are cleared up.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!