ODE45 and 3D arrays

10 views (last 30 days)
Salar
Salar on 27 Jul 2016
Commented: James Tursa on 27 Jul 2016
Hello,
I have a question about 3D-arrays and ODE45. I have a 3D array called "A". "A" contains a set of a matrices at different instants of time, so for example A(:,:,1) will give me a matrix at T_0 . I also have a equation that I need to solve using ODE45. btw all the matrices in A are 6X6. My equation is Phi_dot(:,:,k) = A(:,:,k)*Phi, I'm trying to integrate this, but I keep getting an error on "indices". I would really appreciate it if you could help me with this. Thank you. My code is below :
function [ Phi_dot ] = Phi( A, Phi, t )
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
Phi_dot = A*Phi;
end
Phi_0 = eye(6);
tspan = [0 pi];
opts = odeset('RelTol',1e-12,'AbsTol',1e-12);
for k = 1:i
[t1, Phi] = ode45(@(t,Phi) Phi(A(:,:,k),Phi(:,:,k),t),tspan,Phi_0,opts);
end

Accepted Answer

James Tursa
James Tursa on 27 Jul 2016
Edited: James Tursa on 27 Jul 2016
1) Don't use (:,:,k) subscripting on Phi in your derivative function. ode45 only knows about Phi as a 36-element column vector.
2) Don't use Phi for multiple things! You are using it as the dummy variable in your anonymous function, and as the derivative function name, and as the input variable in your derivative function, and as the return variable from ode45. Use different names for these. Partly because some of these will cause errors, and partly for readability.
3) Make sure your derivative function reshapes the column vector input into a matrix, and then returns a column vector.
E.g.,
[t1, PhiResult] = ode45(@(t,Phi) PhiFun(A(:,:,k),Phi,t),tspan,Phi_0,opts);
:
function [ Phi_dot ] = PhiFun( A, Phi, t )
%UNTITLED2 Summary of this function goes here
% Detailed explanation goes here
Phi_dot = A*reshape(Phi,6,6);
Phi_dot = Phi_dot(:);
end
Then of course you will probably want to add some code to save the desired results from each iteration (the way it is now the results of each iteration overwrite the previous iteration).
  2 Comments
Salar
Salar on 27 Jul 2016
Thank you for your response. My only problem with this is that for my k value, this is taking forever! Would you able to suggest a faster approach maybe I could try? Thank you!
James Tursa
James Tursa on 27 Jul 2016
Do you mean i is large?

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!