solve parallel equations 1

[Amleto]

I am currently using ode45 but am a bit stuck.

What I need to do is basically progress the solution of two instances of the same equation, with different initial conditions, at the same time.

I could just write out the equation again in the function but eventually I may be dealing with 50 of them and I'm sure there must be a better way...

Any help will be most appreciated.

A

 

[MikeyP]

Can you post the odefile for your 2 equation case.

M

--
Dr Michael F Platten

 

[Amleto]

ok, the deal is - I have a few bubbles. They pulsate due to an externally changing applied pressure. The pulsations and the external pressure make the bubbles move, and the following code attempts to model that for two bubbles.

function dr = trans1func(t,r)
global sPeq sPatm sP0 sPv Re We dk epsilon kappa apptype tf a w

% ok, this program considers how two pulsating bubbles translate.
% A lot of code has been duplicated, only with the variable index changed
% and a slight change in code which be removed replacing the change
% with a variable, and re-assigning the variable different values at the
% correct change. so explicitly the equations aren't identical, but
% with a bit of code twidling, the only thing that changes from bubble #1
% to bubble #2 is the indexing.

Cl = 1500 * a * w;

% ------------
% BUBBLE 1
% ------------


% apptype == 0 if the sound field is standing wave, 1 if travelling wave
if apptype == 0  
    sPinf = sP0*(1 + epsilon * cos(dk*r(3) ) * sin(t) );
    sGradPinf = -dk*sP0*epsilon * sin(dk*r(3) ) * sin(t);
    sPinfdot = sP0*epsilon*cos(dk*r(3) ) * cos(t);
elseif apptype == 1 
    sPinf = sP0*(1 + epsilon * sin( dk*r(3) - t ) );
    sGradPinf = dk*sP0*epsilon*cos( dk*r(3) -t );
    sPinfdot = -sP0*epsilon*cos(dk*r(3) - t);
end


% the following calculates the liquid pressure around bubble#1
if r(1)<a
    sPliq = Pgfuncii(r(1),sPeq,kappa) + sPv  ...
           - 4/Re*r(2)./r(1) - 2/We./r(1) ... % adiabatic
        + r(1).*Pgdotfuncii(r(1),r(2),sPeq,kappa)/3./Cgfuncii(r(1),a,1,kappa);
else
    sPliq = Pgfuncii( r(1), sPeq, 1)  + sPv ...
            - 4/Re*r(2)./r(1) - 2/We./r(1) ... % isothermal
            + r(1).*Pgdotfuncii(r(1),r(2),sPeq,1) /3./ Cgfuncii(r(1),a,1,1);
end



dr(1) = r(2); % r(1) = R1,   r(2) = R1dot
% dr(2) = ( 1./r(1) ) * ( sPliq - sPinf - 3.*r(2).^2/2 );

dr(3) = r(4); % r(3) = X1,   r(4) = X1dot
dr(5) = -sGradPinf; % r(5) = Ua,    dr(5) = Ua_dot ~ -GRAD Pinf
dr(4) = (18/Re)*( r(5)-r(4) )./(r(1).^2) + 2*sGradPinf + ...
         3*(r(2)./r(1)).*(r(5)-r(4)) + dr(5) + ...
         2*r(6).^2.*r(7); % this is the line that will
     % only change slightly
     % depending on number of bubbles.
     % the last term is the sum of that product for all bubbles except
     % the one being 'considered' (#1 here)
     

PlmPinf = sPliq - sPinf;
% this term is without the r_tt term
PlmPinfdot = -3*sPeq*kappa*r(2)./r(1).^(3*kappa+1) +...
             r(2)* Pgdotfuncii( r(1),r(2),sP0,kappa )/3/Cgfuncii( r(1),a,1,kappa ) + ...
             r(1)* ( ...
               Pgddotfuncii( r(1),r(2),dr(2),sPeq,kappa )./Cgfuncii( r(1),a,1,kappa ) -...
               Pgdotfuncii( r(1),r(2),sP0,kappa ).* Cgdotfuncii( r(1),r(2),a,sPeq,kappa,1 )/ ...
               Cgfuncii( r(1),a,1,kappa ).^2 ...
             )/3 + (2/We)*r(2)./r(1).^2 - sPinfdot + (4/Re)*( r(2)./r(1) ).^2 - ...
             ( r(5) - r(4) ).*sGradPinf;

 dr(2) = ( PlmPinf + r(1).*PlmPinfdot/Cl - 1.5*r(2).^2 ) ./ ...
         ( r(1) + 4/Cl/Re );

     
     
% ------------
% BUBBLE 2
% ------------     
if apptype == 0
    sPinf = sP0*(1 + epsilon * cos(dk*r(8) ) * sin(t) );
    sGradPinf = -dk*sP0*epsilon * sin(dk*r(8) ) * sin(t);
    sPinfdot = sP0*epsilon*cos(dk*r(8) ) * cos(t);
elseif apptype == 1 
    sPinf = sP0*(1 + epsilon * sin( dk*r(8) - t ) );
    sGradPinf = dk*sP0*epsilon*cos( dk*r(8) -t );
    sPinfdot = -sP0*epsilon*cos(dk*r(8) - t);
end

if r(1)<a
    sPliq = Pgfuncii(r(6),sPeq,kappa) + sPv  ...
           - 4/Re*r(7)./r(6) - 2/We./r(6) ... % adiabatic
        + r(6).*Pgdotfuncii(r(6),r(7),sPeq,kappa)/3./Cgfuncii(r(6),a,1,kappa);
else
    sPliq = Pgfuncii( r(6), sPeq, 1)  + sPv ...
            - 4/Re*r(7)./r(6) - 2/We./r(6) ... % isothermal
            + r(6).*Pgdotfuncii(r(6),r(7),sPeq,1) /3./ Cgfuncii(r(6),a,1,1);
end



dr(6) = r(7); % r(6) = R2,   r(2) = R2dot


dr(8) = r(9); % r(8) = X2,   r(9) = X2dot
dr(5) = -sGradPinf; % r(5) = Ua,    dr(5) = Ua_dot
dr(9) = (18/Re)*( r(5)-r(4) )./(r(1).^2) + 2*sGradPinf + ...
         3*(r(2)./r(1)).*(r(5)-r(4)) + dr(5) + ...
         2*r(1).^2.*r(2); % this is the line that will
     % only change slightly
     % depending on number of bubbles.
     % the last term is the sum of that product for all bubbles except
     % the one being 'considered' (#2 here)

PlmPinf = sPliq - sPinf;
% this term is without the r_tt term
PlmPinfdot = -3*sPeq*kappa*r(7)./r(6).^(3*kappa+1) +...
             r(7)* Pgdotfuncii( r(6),r(7),sP0,kappa )/3/Cgfuncii( r(6),a,1,kappa ) + ...
             r(7)* ( ...
               Pgddotfuncii( r(6),r(7),dr(7),sPeq,kappa )./Cgfuncii( r(6),a,1,kappa ) -...
               Pgdotfuncii( r(6),r(7),sP0,kappa ).* Cgdotfuncii( r(6),r(7),a,sPeq,kappa,1 )/ ...
               Cgfuncii( r(6),a,1,kappa ).^2 ...
             )/3 + (2/We)*r(7)./r(6).^2 - sPinfdot + (4/Re)*( r(7)./r(6) ).^2 - ...
             ( r(10) - r(9) ).*sGradPinf;

 dr(7) = ( PlmPinf + r(6).*PlmPinfdot/Cl - 1.5*r(7).^2 ) ./ ...
         ( r(6) + 4/Cl/Re );

% -------------
dr=dr';

Hope it's well enough annotated

 

[MikeyP]

Just noticed that there is no definition for dr(10) in bubble 2. I presume that the dr(5) line in bubble 2 should read dr(10).

Well, the obvious - but probably not the most efficient - thing to do is to put everything in a for loop and build up a vector of "dr's" by (as you suggest) changing the indexing relative to the loop variable. However this will be very slow for large numbers of bubbles.

A better way is to vectorise the code with respect to r (you don't need to vectorise with respect to t if you are using ode45). So for 3 bubbles say, r will be a 15 element vector with 1, 6 and 11 being the first r variable for each bubble and 2, 7, 12 bing the second etc. Each dr will then be of the form

bubbles = 3; % number of bubbles
dr(1:5:5*(bubbles-1)+1) = some_function(r(1:5:5*(bubbles-1)+1));
dr(2:5:5*(bubbles-1)+2) = another_function(r(2:5:5*(bubbles-1)+2);
% etc.

For example, your 3rd variable in the derivative will change from

dr(3) = r(4);

to

dr(3:5:5*(bubbles-1)+3) = r(4:5:5*(bubbles-1)+3);

Of course this requires that your various functions *func* and *funcdot* are also vectorised with respect to r.

The 4th variable in the derivative is a little trickier (but not much) because of the last term. Sum the last term for all bubbles and then subtract the term for the "current" bubble) thus:

sum_prod = sum(2.*r(1:5:5*(bubbles-1)+1).^2.*r(2:5:5*(bubbles-1)+2)); % the sum of product for ALL bubbles (a scalar)

prod_bubble = (2.*r(1:5:5*(bubbles-1)+1).^2.*r(2:5:5*(bubbles-1)+2); % the product for each individual bubble (a vector)

dr(4:5:(bubbles-1)+4) = (18/Re ... blah blah ... + (sum_prod - prod_bubble);

M







--
Dr Michael F Platten

 

[Amleto]

hi, thanks a lot for the help.

dr(5) is a universal quantity and doesn't depend on the bubbles so maybe I should make this dr(1), and then the start of the bubble quantities would be 2,8,12 etc?

Would this require all r(x:5:z) to be changed to r(x:4:z)?


I would have thought

dr(3:5:5*(bubbles-1)+3) = r(4:5:5*(bubbles-1)+[COLOR=red]3[/color]);

should be

dr(3:5:5*(bubbles-1)+3) = r(4:5:5*(bubbles-1)+[COLOR=red]4[/color]);

?

thanks again,

Tom