Using a Subroutine to Create a Plot

I have the following subroutine which takes a point with radius 1 and then finds a velocity vector at that point.

p = [0.282 0 0.959]';          %specify point where velocity is to be evaluated%

[v]=StokesletVelocity(N,sing,p,Kn,f);

 function[v]=StokesletVelocity(N,sing,p,Kn,f)
 for i=1:N
  rij= (p - sing(1:3,i));
  FF(1:3, 3*i-2:3*i) =((1/(8*pi*Kn))*((eye(3)/norm(rij)) +  (rij*rij')/norm(rij)^3))';          %subroutine to find velocity at the point %
 end
 
  v=FF*f;
end

What I would like to do is plug in points to the subroutine with a range of radii going from 1 to 10 and get velocities for all of them. I was thinking of doing this in steps of 0.2 at a time, so the next p to plug into the subroutine would be

[0.2*0.282 0 0.2*0.959]

the next would be

 [0.4*0.282 0 0.4*0.959]

and so on up until a radius of 10. This will produce an array of velocity vectors and I then need to find the radial component for each one and plot the radial component of each velocity vector against the corresponding value for the radius ie. the number which is mutliplying the numbers in the position vector for the point. Is it possible to use the subroutine I have to generate the array of velocity vectors and then plot radial component against radius in this way? Let me know if it is not clear what I mean.

Answer 1

[Jon]

I think that what you describe could be implemented with a simple loop. So make a vector of radii that you want to evaluate, set up a for loop that goes through the radii values, calculates p, calls your stokes function, and saves the results. Once you exit the loop you extract your velocity component and plot it agains the radius.

What are the dimensions of the variable v, that is returned by StokesletVelocity?

For illustrative purposes I will assume it is 3 by 1.

Then you could do something like this

% define the radius values
radius = 0.2:0.2:10 % radius goes from 0.2 to 10 in increments of 0.2

% define initial position
p = [0.282 0 0.959]

% you will need to assign numerical values for these variables
% I put them outside the loop assuming that they are not functions of the radius
% you will get an error on these lines until you fill in something on the
% right hand side of these equations.
N = 
sing = 
Kn = 
f = 

% preallocate array to hold result
numRadius = length(radius); % number of radius values
V = zeros(3,numRadius); % columns are results for each radius value, rows are velocity components
% loop through the values of radius calculating velocities
for k = 1:numRadius
    p = radius(k)*p;
    
    % evaluate the velocity and store it
    V(:,k) = StokesletVelocity(N,sing,p,Kn,f);
    
end

% plot one of the components versus radius
plot(radius,V(1,:))

You may have to tweak the above slightly, as there are aspects of your problem that I couldn't fully understand from your description. Hopefully this gets you close enough that you can adapt it from here.

[Tom]

I have used something like this to get the array of velocity vectors at the points along the radial line, I now need the radial and polar component of each velocity to be plotted against the corresponding radius, I suppose this might be possible with the sph2cart type functions in MATLAB?

So perhaps I did not explain this properly, but in the array of velocity vectors each velocity is in Cartesian components, and I then need to consider each velocity, extract the radial component and plot against radius, then do the same for the polar component of each velocity vector in the array.

[Jon]

There is also cart2sph to go from cartesian to spherical. I'm not quite clear whether you have spherical coordinates, or some other polar representation. Anyhow are you clear about how to add this to your code now?

[Tom]

Hi Jon, yes I think I confused myself a bit, the initial point

p = [0.282 0 0.959]

is in Cartesian coordinates so the method of multiplying directly by the radius might not work, although I can convert it to spherical polar and start from there.

I think there is something wrong with the loop you have suggested, as the graph I am creating starts at 1 as I expect but then trails off to quickly, I think this is because the loop takes the base radius and multiplies by 1.2, then in the next loop it multiplies by 1.4 x 1.2, then in the next by 1.6 x 1.4 x 1.2, is there a way of fixing this?