ode45 code issue solving for polar coordinates

Okay so I trying to get ode45 to work. To do this I established a function V which represents the 3 velocities of a fluid. I establish all the intial conditions in the function as well as the varialbe. I want the ode45 function to output the 3 variables (r,T,z) Eventually I am going to use pol2cart to convert and plot the function in the x,y,z plan. Can anyone assist in pointing out were I am going wrong.

b = 2;
aspect = 0.5;
c = 3;
a = 2;
r = 2;
T = 0.1;
z = 5;
t = 0;
n = 50000;
Init_cond = 1;
tspan = [0 n];

% The values above are not the final values I am currently reseaching 
% What they represent and trying to determine an appropate model
% They will be updated soon


[r,T,z] = ode45(@odeFun, tspan, Init_cond);
[x,y,z] = pol2cart(r,T,z);
plot3(x,y,z);
% using the matlab funciton ode45 I want it to output the variables r,T,z
% input function is odeFun which is equal to the velocity functions
% origianlly provided
% Initial condition is just so for all variables so r, T and z = 0 when
% time equals zero


function V = odeFun(r,T,z)
b = 2;
aspect = 0.5;
c = 3;
a = 2;   
r(1) = 0; % initial conditions
T(1) = 0; % initial conditions
z(1) = 0; % initial conditions
V = zeros(3,2,1);
V(1) = @(r,T,z) r*(z-b)/aspect; % Velocity in the radius
V(2) = @(r,T,z) c; % Velocity in the Theta
V(3) = @(r,T,z) z*(1-a*r); % Velocity in the z direction
end

[Philip Binaco]

I made a couple small changes that are not working either but I am my thought is that I was originally trying to skip a step would i need to do 3 different ode45 to find the variables [t,r] [t,T] and [t,z]. Then I could potentially do a pol2cart of and plot the the [x,y,z]

b = 2;
aspect = 0.5;
c = 3;
a = 2;
r = 2;
T = 0.1;
z = 5;
t = 0;
n = 50000;
Init_cond = 1;
tspan = [0 n];

% The values above are not the final values I am currently reseaching 
% What they represent and trying to determine an appropate model
% They will be updated soon


[t,r] = ode45(@odeFun, tspan, Init_cond);
[t,T] = ode45(@odeFun, tspan, Init_cond);
[t,z] = ode45(@odeFun, tspan, Init_cond);
% using the matlab funciton ode45 I want it to output the variables r,T,z
% input function is odeFun which is equal to the velocity functions
% origianlly provided
% Initial condition is just so for all variables so r, T and z = 0 when
% time equals zero



function dVdt = odeFun(r,T,z);
b = 2;
aspect = 0.5;
c = 3;
a = 2;   
r(1) = 0; % initial conditions
T(1) = 0; % initial conditions
z(1) = 0; % initial conditions
dVdt = zeros(3,2,1);
dVdt(1) = @(r,T,z) r*(z-b)/aspect; % Velocity in the radius
dVdt(2) = @(r,T,z) c; % Velocity in the Theta
dVdt(3) = @(r,T,z) z*(1-a*r); % Velocity in the z direction
end

Answer 1

[Philip Binaco]

Finally got a working code here

clear all

t = 0;
n = 1000;
Init_cond = zeros(3,1);
Init_cond(2) = 0.4; % initial condition for r(1) = 0.1
Init_cond(1) = 4; % initial condition for T(1) = 0.3
Init_cond(3) = 100; % inital condition for z(1) = 100
tspan = [t n];

% The values above are not the final values I am currently reseaching 
% What they represent and trying to determine an appropate model
% They will be updated soon


[t,Polar] = ode45(@odeFun, tspan, Init_cond);
[cart(:,1),cart(:,2),cart(:,3)] = pol2cart(Polar(:,1),Polar(:,2),Polar(:,3));

%plot3(Polar(:,1), Polar(:,2), Polar(:,3))
plot3(cart(:,1), cart(:,2), cart(:,3))
% subplot(2,1,1)
% plot(Polar(:,2),Polar(:,3))
% xlabel('r(t)')
% ylabel('z(t)')

% subplot(2,1,2)
% polarplot(Polar(:,1),Polar(:,2),'r-');

% subplot(3,1,3);
% plot3(cart(:,1),cart(:,2),cart(:,3),'r-');
% xlabel('x(t)')
% ylabel('y(t)')
% zlabel('z(t)')

% using the matlab funciton ode45 I want it to output the variables r,T,z
% input function is odeFun which is equal to the velocity functions
% origianlly provided
% Initial condition is just so for all variables so r, T and z = 0 when
% time equals zero



function dVdt = odeFun(t,y)
b =20;
aspect = 60;
c = 5;
a = 5;
r = y(2);
T = y(1);
z = y(3);
dVdt = zeros(3,1);
dVdt(2) = r*(z-b)/aspect; % Velocity in the radius  
dVdt(1) = c; % Velocity in the Theta
dVdt(3) = z*(1-a*r); % Velocity in the z direction
end

Answer 2

[Walter Roberson]

dVdt(3) = @(r,T,z) z*(1-a*r); % Velocity in the z direction

That attempts to create an anonymous function and store the function handle as the third element of dVdt . However, it is not permitted to create ()-indexed arrays of function handles, so if you needed to create an array of function handles you would need

dVdt = cell(3,1);
dVdt{1} = @(r,T,z) r*(z-b)/aspect; % Velocity in the radius
dVdt{2} = @(r,T,z) c; % Velocity in the Theta
dVdt{3} = @(r,T,z) z*(1-a*r); % Velocity in the z direction

... however, what you return from the ode function must be strictly datatype single() or datatype double() .

What you return from the ode function must be a column vector that has the same number of elements as your initial condition. Your initial condition is a scalar, so unless you change your initial conditions, you would need to return a scalar.

[Walter Roberson]

I want the ode45 function to output the 3 variables (r,T,z)

Your approach is wrong. You need to bundle the three different variables into a single vector, and pass in a vector of three initial conditions, and your ode function needs to return a vector of derivatives the same size as your initial conditions.

r = 2;
T = 0.1;
z = 5;

Are those your initial conditions for r, T, z?

r(1) = 0; % initial conditions
T(1) = 0; % initial conditions
z(1) = 0; % initial conditions

or are those your initial conditions for r, T, z ?

[Philip Binaco]

Hey Walter,

Thank you for your response I updated the code slightly but:

r(1) = 2;
T(1) = 0.1;
z(1) = 5;

^^^^^^^Are supposed to be the initial condition. But what I did after update my code was create a Init_cond = cell(3,1) this is were i am storing my r(1) = 2, T(1) = 0.1 and z(1) = 5

Attached below are my

t = 0;
n = 50000;
Init_cond = cell(1,3);
Init_cond{1} = 2; % initial condition for r(1) = 2
Init_cond{2} = 0.1; % initial condition for T(1) = 0.1
Init_cond{3} = 5; % inital condition for z(1) = 5
tspan = [0 n];

% The values above are not the final values I am currently reseaching 
% What they represent and trying to determine an appropate model
% They will be updated soon


[r,T,z] = ode45(@odeFun, tspan, Init_cond);

% using the matlab funciton ode45 I want it to output the variables r,T,z
% input function is odeFun which is equal to the velocity functions
% origianlly provided
% Initial condition is just so for all variables so r, T and z = 0 when
% time equals zero



function dVdt = odeFun(r,T,z);
    b = 2;
    aspect = 0.5;
    c = 3;
    a = 2;
    dVdt = cell(3,1);
    dVdt{1} = @(r,T,z) r*(z-b)/aspect; % Velocity in the radius
    dVdt{2} = @(r,T,z) c; % Velocity in the Theta
    dVdt{3} = @(r,T,z) z*(1-a*r); % Velocity in the z direction
    %dVdt = zeros(3,2,1);
    %dVdt(1) = @(r,T,z) r*(z-b)/aspect; % Velocity in the radius
    %dVdt(2) = @(r,T,z) c; % Velocity in the Theta
    %dVdt(3) = @(r,T,z) z*(1-a*r); % Velocity in the z direction
end

[Walter Roberson]

As I wrote above,

... however, what you return from the ode function must be strictly datatype single() or datatype double() .

Returning a cell array of anonymous functions is not permitted .

The ode*() routines do numeric integration. They figure out a step size and they propose a series of boundary conditions based upon continuity models, and they test to see whether the proposed positions give values within tolerance for the model. If the proposal is accepted, then they take the returned numeric values (which are instantaneous derivatives), multiply them by the step size, and add those to the previous boundary conditions to get new boundary conditions.

The routines do not expect anonymous functions to be returned, and will not call any returned anonymous functions to determine numeric values: the routines test whether isfloat() the returned values, and if not then error() .