Summation of Function Handle Cells

Hello.

I am trying to sum up each individual cell that is produced by a function handle in a for loop that runs from.

My handle looks like this:

for a = 1 : 11;
    f{a} = @(u)1
    for b = setdiff(1:11, a); 
        f{a} = @(u)1.*G.*mass_a.*mass_b.*r./(norm_r^3)
        %need to sum up all 11 cells produced here
    end
end

The function handle spits out something like this (I have only put 2 cells, there are 11):

{@(u)1.*G.*mass_a.*mass_b.*r./(norm_r^3)}    {@(u)1.*G.*mass_a.*mass_b.*r./(norm_r^3)}

I need to sum all of the cells from above up.

It seems like the sum function on matlab does not work. I have seen a couple of instances of loops being used, but can't seem to be able to apply them correctly.

Thanks

[Geoff Hayes]

AeroLad - I don't understand why in the outer loop you do

f{a} = @(u)1

and in the inner loop you overwrite this function handle with

f{a} = @(u)1.*G.*mass_a.*mass_b.*r./(norm_r^3)

? Is this intentional? f is a cell array of function handles. At what point do you want to evaluate the function for a some value? Please clarify.

[AeroLad]

Hi Geoff, this is intentional. To give a bit of background, I am writing code to simulate the n-body problem.

%pp is a 11x7 struct
%Sun mass = pp(1:1, 1:1) Sun position = pp(1:1, 2:4)  Sun Velocity = pp(1:1, 5:7)
%Mercury mass = pp(2:2, 1:1) Sun position = pp(2:2, 2:4)  Sun Velocity = pp(2:2, 5:7)
%Similarly for rest of planets

[XX,No_Planets] = size(Planets)
%size = (1,11)


t0 = 0; %intital time (days)
T = 365*0.5; %final time
dt = 10; %timestep (days) 
for t = T/dt

    for a = 1 : No_Planets;
      
        disp('Planet a:'); disp(a)
        mass_a = pp(a:a,1:1);
        r_a = [pp(a:a,2:2),pp(a:a,3:3),pp(a:a,4:4)];
        v_a = [pp(a:a,5:5),pp(a:a,6:6),pp(a:a,7:7)];
        
        f{a} = @(u)1
        
   %Loop B compares all of the planets/ bodies (except a)with planet/body a

        for b = setdiff(1:No_Planets, a); 
            
            disp('Planet b:'); disp(b)
            %Positions of a planet that is not the one chosen in 'a'
            mass_b = pp(b:b,1:1);
            r_b = [pp(b:b,2:2),pp(b:b,3:3),pp(b:b,4:4)]
            v_b = [pp(b:b,5:5),pp(b:b,6:6),pp(b:b,7:7)]
            

            r =  r_b -r_a;
            r = r'
            norm_r = norm(r);
            
            f{a} = @(u)1.*G.*mass_a.*mass_b.*r./(norm_r^3)
        
            z = zeros(1,11)  %% Attempt to sum up the function handle. Probably very wrong
            for k = 1:b
                fh = @(u) fh + f(k,:)
            end
            
        end
    end
%fv =arrayfun( @(u) fk(u), [1;11] )


    end

The aim is to get equations of motion described by .

The vectors r_a, r_b, v_a, v_b produce something like this: 1.0e+09 * [2.865 3.383 0.1924].

In the code, r =

The summation loop will then for each row basically hold 'a' (i in the formula) constant, while taking all values of b (j in the formula), then move on to the next a, and take all 10 other values of b etc...

The outermost loop (T/dt) is still a work in progress

Thanks

[Walter Roberson]

f{a} = @(u)1

It is usually better to use

f{a} = @(u) ones(size(u))

and likewise for equations such as

f{a} = @(u)1.*G.*mass_a.*mass_b.*r./(norm_r^3)

it is best to multiply by ones(size(u))

It might possibly not matter for your situation, but it is best to get into the habit of this.

The reason it makes a difference is that if u is a non-scalar then @(u)1 returns a scalar rather than being the same size as u, and that is a problem for vectorized use in functions such as integral().

[AeroLad]

Thanks. Incorporated. Would you know how I could sum the cells of the function handle? Spent hours and hours looking, but not quite found the right way. Maybe I am just looking wrong. I am by no means a matlab expert.

[James Tursa]

Why all the function handles? Why not just calculate the forces directly inside your loops and add them up? Simpler and quicker I would think.

[AeroLad]

I'm probably being dumb here, but when I do this:

f{a} = @(u)ones(size(u))

I get the following error:

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

[AeroLad]

@James, how exactly does on go about doing that? A bit lost on the process. It would naturally fall into the b loop, but how do I iteratively add all the terms of b?

I am trying to get something like this when a= 1:

f = @(u) G*m2*(r_b (iteration 1)-r_a )/r^3 + G*m2*(r_b (iteration 2)-r_a )/r^3 + G*m2*(r_b (iteration 3)-r_a )/r^3...

Any efficient method of doing so would be much appreciated

[Walter Roberson]

You have an existing f variable in your workspace that is not a cell array.

... Using functions helps reduce this kind of problem.

[AeroLad]

Nail on the head! Thanks

[James Tursa]

@AeroLad: My point is that the way you are currently building your function handles the radius values r_a and r_b are only snapshots of the current r_a and r_b for that iteration. I.e., all of the function handles you are building and summing are only good for one iteration of your calculations. E.g., If you are using this in an iterative integration scheme, you would be forced to completely rebuild all of your function handles from scratch at each step in the iteration process. Is that what you intend? It seems to me that building these function handles doesn't buy you anything. Just calculating the forces directly and summing them would be quicker. But I don't know the overall goal you have.

[AeroLad]

Ah. I think I see what you mean. I will look into calculating them directly. It will also most certainly make things faster. Thanks

Answer 1

[Walter Roberson]

fsum = @(x) sum(cellfun(@(F) F(x), f));

[AeroLad]

Many thanks! All the best and stay safe!

[Omid Feizi]

It was my question too, your code help me a lot.