Advanced Search

Doug Meade

3088 Reputation

23 Badges

20 years, 336 days
Contact Doug Meade

 

Doug

---------------------------------------------------------------------
Douglas B. Meade  <><
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu

MaplePrimes Activity


These are replies submitted by Doug Meade

jerks don't jump...

Doug Meade 3088
February 20 2008
0 0

I thought the code, and graphs, in your worksheet looked pretty familiar. I'm glad you are finding them useful.

The way I addressed this issue was to model the deployment over a short time interval. To keep the jerk continuous I used a cubic function. In the end, this has very little affect on the final landing time or velocity - it just gets you through the deployment without being torn apart.

I believe I did solve the problem of finding the latest time to deploy the chute while keeping the impact velocity below a certain threshold. I believe I did this iteratively. You might be able to do something with a shooting method.

Because the ODEs are autonomous and the velocity at deployment is essentially terminal velocity, I believe you can work out the minimal time needed to bring the velocity from (pre-deployment) terminal velocity your threshold, then work backwards from here. This will get you close; a couple iterations should lead to a good estimate.

That was a fun project. I still like to talk about it.

Let me know if I can be of further help.

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

jerks don't jump...

Doug Meade 3088
February 20 2008

I thought the code, and graphs, in your worksheet looked pretty familiar. I'm glad you are finding them useful.

The way I addressed this issue was to model the deployment over a short time interval. To keep the jerk continuous I used a cubic function. In the end, this has very little affect on the final landing time or velocity - it just gets you through the deployment without being torn apart.

I believe I did solve the problem of finding the latest time to deploy the chute while keeping the impact velocity below a certain threshold. I believe I did this iteratively. You might be able to do something with a shooting method.

Because the ODEs are autonomous and the velocity at deployment is essentially terminal velocity, I believe you can work out the minimal time needed to bring the velocity from (pre-deployment) terminal velocity your threshold, then work backwards from here. This will get you close; a couple iterations should lead to a good estimate.

That was a fun project. I still like to talk about it.

Let me know if I can be of further help.

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

parachute...

Doug Meade 3088
February 20 2008
0 0

Look at the graph in your worksheet. The velocity has been (essentially) constant for more than 3 minutes. Unless you deploy the parachute so late that impact occurs before terminal velocity is reached, the impact velocity is essentially independent of the deployment time.

But first, are you sure the skydiver is not torn to pieces long before impact?

Compute the force on the skydiver at the instant of deployment. This can be done without solving a differential equation, just compute the force before and after deployment (position and velocity are continuous here, but not acceleration). I find the jerk is 6G. Not as bad as most textbook problems from a decade ago.

I did some modeling of this problem about 10 years ago. Here are two entries from my list of publications:

Back to your original question: This is a difficult problem for Maple to solve. You really need an exact solution (which is not difficult to write) and then some messy manipulations. Before I invest too much time in this, why do you want to know this relationship? Maybe there is another way to get this from the problem.

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

parachute...

Doug Meade 3088
February 20 2008

Look at the graph in your worksheet. The velocity has been (essentially) constant for more than 3 minutes. Unless you deploy the parachute so late that impact occurs before terminal velocity is reached, the impact velocity is essentially independent of the deployment time.

But first, are you sure the skydiver is not torn to pieces long before impact?

Compute the force on the skydiver at the instant of deployment. This can be done without solving a differential equation, just compute the force before and after deployment (position and velocity are continuous here, but not acceleration). I find the jerk is 6G. Not as bad as most textbook problems from a decade ago.

I did some modeling of this problem about 10 years ago. Here are two entries from my list of publications:

Back to your original question: This is a difficult problem for Maple to solve. You really need an exact solution (which is not difficult to write) and then some messy manipulations. Before I invest too much time in this, why do you want to know this relationship? Maybe there is another way to get this from the problem.

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Let's think about it first...

Doug Meade 3088
February 20 2008
0 0

This is a messy integral. Even though the interval is reasonable, [0.35,1], the Bessel function is being evaluated on an interval that it 50001 times larger, [17500.35,50001]. Over this interval, the integrand oscillates quite a few times. An exact value for this integral is unlikely and the oscillations make it difficult to compute with any accuracy.

I am not sure why you were trying to plot the value of a definite integral. Here are a few steps to get a better view on this problem.

F := Int(BesselJ(0, 50001*x)*x*exp(I*(355*x^2*1/2)), x = .35 .. 1);
f := op(1,F);
fR := evalc( Re(f) );
plot( fR, x=0.35..1 );
evalf( Int( fR, x=0.35..1 ) );  # returns unevaluated

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Let's think about it first...

Doug Meade 3088
February 20 2008

This is a messy integral. Even though the interval is reasonable, [0.35,1], the Bessel function is being evaluated on an interval that it 50001 times larger, [17500.35,50001]. Over this interval, the integrand oscillates quite a few times. An exact value for this integral is unlikely and the oscillations make it difficult to compute with any accuracy.

I am not sure why you were trying to plot the value of a definite integral. Here are a few steps to get a better view on this problem.

F := Int(BesselJ(0, 50001*x)*x*exp(I*(355*x^2*1/2)), x = .35 .. 1);
f := op(1,F);
fR := evalc( Re(f) );
plot( fR, x=0.35..1 );
evalf( Int( fR, x=0.35..1 ) );  # returns unevaluated

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Now I understand (better)...

Doug Meade 3088
February 18 2008
0 0

OK, I think I understand your question, and see why it takes so long.

Use the suggestions already posted to get an explicit formula for ETa, then replace Pna with Pa-Inta. This is going to be very messy. This can now be evaluated for specific values of Pa. A plot is going to take quite a long time to create due to the number of improper integrals that need to be evaluated numerically.

restart;
Inta := Pa*Int(Pm*(1-exp(-1.220/(Pm^.45-.2*Pm^.69+0.2e-1*Pm)))*exp(-Pm/(0.97e-1*Pa+8.57))/(0.97e-1*Pa+8.57), Pm = 0 .. infinity)/(0.97e-1*Pa+8.57):
ETa:= 9.285051068*Int(9.285051068*min(.5676676416*Pnm, 50.)*exp(-9.285051068*Pnm/Pna)/Pna, Pnm = 0. .. Float(infinity)):
ETa2 := value(convert( ETa, piecewise, Pnm )) assuming Pna>0:
ETa3 := eval( ETa2, Pna=Pa-Inta ):
ETa4 := unapply( evalf(ETa3), Pa ):
evalf( ETa4( 0.1 ) );
                                0.03277703754

This calculation took about 1 CPU minute

My suggestion is to find a good (enough) approximation to the improper integrals - truncate the interval and/or approximate the integrand.

Good luck!

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Now I understand (better)...

Doug Meade 3088
February 18 2008

OK, I think I understand your question, and see why it takes so long.

Use the suggestions already posted to get an explicit formula for ETa, then replace Pna with Pa-Inta. This is going to be very messy. This can now be evaluated for specific values of Pa. A plot is going to take quite a long time to create due to the number of improper integrals that need to be evaluated numerically.

restart;
Inta := Pa*Int(Pm*(1-exp(-1.220/(Pm^.45-.2*Pm^.69+0.2e-1*Pm)))*exp(-Pm/(0.97e-1*Pa+8.57))/(0.97e-1*Pa+8.57), Pm = 0 .. infinity)/(0.97e-1*Pa+8.57):
ETa:= 9.285051068*Int(9.285051068*min(.5676676416*Pnm, 50.)*exp(-9.285051068*Pnm/Pna)/Pna, Pnm = 0. .. Float(infinity)):
ETa2 := value(convert( ETa, piecewise, Pnm )) assuming Pna>0:
ETa3 := eval( ETa2, Pna=Pa-Inta ):
ETa4 := unapply( evalf(ETa3), Pa ):
evalf( ETa4( 0.1 ) );
                                0.03277703754

This calculation took about 1 CPU minute

My suggestion is to find a good (enough) approximation to the improper integrals - truncate the interval and/or approximate the integrand.

Good luck!

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

What did you expect?...

Doug Meade 3088
February 16 2008
0 0

It is not possible to solve the first equations for beta, in closed form. But, if you take the final result and expand it as a series you find:

eq1:= tan(beta)= -(2/alpha)*beta;
                                          2 beta
                            tan(beta) = - ------
                                          alpha 
eq2:= beta=sqrt(lambda*alpha - (alpha/2)^2);
                                                     (1/2)
                          1 /                      2\     
                   beta = - \4 lambda alpha - alpha /     
                          2                               
eq3:=eval(eq2,beta = solve(eq1,beta));
                                                                 (1/2)
                                      1 /                      2\     
       RootOf(2 _Z + tan(_Z) alpha) = - \4 lambda alpha - alpha /     
                                      2                               
eq4:=solve(eq3,lambda);
              1       /                                     2\
              - alpha \1 + tan(RootOf(2 _Z + tan(_Z) alpha)) /
              4                                               
series( eq4, alpha=0 );
                            1          /     13\
                            - alpha + O\alpha  /
                            4

This is the series with the default order: 6. Increasing the order to 20 and 200 gives:

series( eq4, alpha=0, 20 );
                            1          /     41\
                            - alpha + O\alpha  /
                            4                   
series( eq4, alpha=0, 200 );
                            1          /     401\
                            - alpha + O\alpha   /
                            4

Based on this, it certainly does appear as though you should have lambda=alpha/4.

The problem is that RootOf is always satisfied by _Z=0, but this is not the only solution. You need to get Maple to return a non-zero solution to the 2*_Z+tan(_Z)*alpha=0. To get a handle on approximate solutions, you might want to look at a series expansion for this equation (in _Z) and use these solutions as approximations to the non-zero solutions to 2*_Z+tan(_Z)*alpha=0.

Here is how this might proceed:

eq5 := solve( eq2, lambda );
                                   2         2
                              alpha  + 4 beta 
                              ----------------
                                  4 alpha     
eq6 := eval( eq5, beta = solve(eq1,beta) );
                       2                                 2
                  alpha  + 4 RootOf(2 _Z + tan(_Z) alpha) 
                  ----------------------------------------
                                  4 alpha                 
series( eq6, alpha );
                            1          /     11\
                            - alpha + O\alpha  /
                            4                   
eq7 := series( 2*_Z+tan(_Z)*alpha, _Z=0 );
                             1         3   2          5    /  6\
            (2 + alpha) _Z + - alpha _Z  + -- alpha _Z  + O\_Z /
                             3             15                   
sol := solve( eval(eq7,O=0)=0, _Z );
                                                       (1/2)  
            /                                    (1/2)\       
            |            /         2            \     |       
          1 |  5 alpha - \-95 alpha  - 240 alpha/     |       
          - |- ---------------------------------------|     , 
          2 \                   alpha                 /       

                                                           (1/2)  
                /                                    (1/2)\       
                |            /         2            \     |       
              1 |  5 alpha - \-95 alpha  - 240 alpha/     |       
            - - |- ---------------------------------------|     , 
              2 \                   alpha                 /       

                                                         (1/2)  
              /                                    (1/2)\       
              |            /         2            \     |       
            1 |  5 alpha + \-95 alpha  - 240 alpha/     |       
            - |- ---------------------------------------|     , 
            2 \                   alpha                 /       

                                                           (1/2)   
                /                                    (1/2)\        
                |            /         2            \     |        
              1 |  5 alpha + \-95 alpha  - 240 alpha/     |        
            - - |- ---------------------------------------|     , 0
              2 \                   alpha                 /        
eq8 := seq( simplify(eval( eq5, beta=B )), B=[sol] );
          3              (1/2)                         (1/2)  
     alpha  - 5 alpha + 5      (-alpha (19 alpha + 48))       
     -------------------------------------------------------, 
                                   2                          
                            4 alpha                           

            3              (1/2)                         (1/2)  
       alpha  - 5 alpha + 5      (-alpha (19 alpha + 48))       
       -------------------------------------------------------, 
                                     2                          
                              4 alpha                           

            3              (1/2)                         (1/2)  
       alpha  - 5 alpha - 5      (-alpha (19 alpha + 48))       
       -------------------------------------------------------, 
                                     2                          
                              4 alpha                           

            3              (1/2)                         (1/2)         
       alpha  - 5 alpha - 5      (-alpha (19 alpha + 48))       1      
       -------------------------------------------------------, - alpha
                                     2                          4      
                              4 alpha

You will have to decide which, if any, of these approximations is relevant for your problem.

plot( [eq8], alpha=-2*Pi..2*Pi, view=[DEFAULT,-10..10] );

Hoping this is helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

What did you expect?...

Doug Meade 3088
February 16 2008

It is not possible to solve the first equations for beta, in closed form. But, if you take the final result and expand it as a series you find:

eq1:= tan(beta)= -(2/alpha)*beta;
                                          2 beta
                            tan(beta) = - ------
                                          alpha 
eq2:= beta=sqrt(lambda*alpha - (alpha/2)^2);
                                                     (1/2)
                          1 /                      2\     
                   beta = - \4 lambda alpha - alpha /     
                          2                               
eq3:=eval(eq2,beta = solve(eq1,beta));
                                                                 (1/2)
                                      1 /                      2\     
       RootOf(2 _Z + tan(_Z) alpha) = - \4 lambda alpha - alpha /     
                                      2                               
eq4:=solve(eq3,lambda);
              1       /                                     2\
              - alpha \1 + tan(RootOf(2 _Z + tan(_Z) alpha)) /
              4                                               
series( eq4, alpha=0 );
                            1          /     13\
                            - alpha + O\alpha  /
                            4

This is the series with the default order: 6. Increasing the order to 20 and 200 gives:

series( eq4, alpha=0, 20 );
                            1          /     41\
                            - alpha + O\alpha  /
                            4                   
series( eq4, alpha=0, 200 );
                            1          /     401\
                            - alpha + O\alpha   /
                            4

Based on this, it certainly does appear as though you should have lambda=alpha/4.

The problem is that RootOf is always satisfied by _Z=0, but this is not the only solution. You need to get Maple to return a non-zero solution to the 2*_Z+tan(_Z)*alpha=0. To get a handle on approximate solutions, you might want to look at a series expansion for this equation (in _Z) and use these solutions as approximations to the non-zero solutions to 2*_Z+tan(_Z)*alpha=0.

Here is how this might proceed:

eq5 := solve( eq2, lambda );
                                   2         2
                              alpha  + 4 beta 
                              ----------------
                                  4 alpha     
eq6 := eval( eq5, beta = solve(eq1,beta) );
                       2                                 2
                  alpha  + 4 RootOf(2 _Z + tan(_Z) alpha) 
                  ----------------------------------------
                                  4 alpha                 
series( eq6, alpha );
                            1          /     11\
                            - alpha + O\alpha  /
                            4                   
eq7 := series( 2*_Z+tan(_Z)*alpha, _Z=0 );
                             1         3   2          5    /  6\
            (2 + alpha) _Z + - alpha _Z  + -- alpha _Z  + O\_Z /
                             3             15                   
sol := solve( eval(eq7,O=0)=0, _Z );
                                                       (1/2)  
            /                                    (1/2)\       
            |            /         2            \     |       
          1 |  5 alpha - \-95 alpha  - 240 alpha/     |       
          - |- ---------------------------------------|     , 
          2 \                   alpha                 /       

                                                           (1/2)  
                /                                    (1/2)\       
                |            /         2            \     |       
              1 |  5 alpha - \-95 alpha  - 240 alpha/     |       
            - - |- ---------------------------------------|     , 
              2 \                   alpha                 /       

                                                         (1/2)  
              /                                    (1/2)\       
              |            /         2            \     |       
            1 |  5 alpha + \-95 alpha  - 240 alpha/     |       
            - |- ---------------------------------------|     , 
            2 \                   alpha                 /       

                                                           (1/2)   
                /                                    (1/2)\        
                |            /         2            \     |        
              1 |  5 alpha + \-95 alpha  - 240 alpha/     |        
            - - |- ---------------------------------------|     , 0
              2 \                   alpha                 /        
eq8 := seq( simplify(eval( eq5, beta=B )), B=[sol] );
          3              (1/2)                         (1/2)  
     alpha  - 5 alpha + 5      (-alpha (19 alpha + 48))       
     -------------------------------------------------------, 
                                   2                          
                            4 alpha                           

            3              (1/2)                         (1/2)  
       alpha  - 5 alpha + 5      (-alpha (19 alpha + 48))       
       -------------------------------------------------------, 
                                     2                          
                              4 alpha                           

            3              (1/2)                         (1/2)  
       alpha  - 5 alpha - 5      (-alpha (19 alpha + 48))       
       -------------------------------------------------------, 
                                     2                          
                              4 alpha                           

            3              (1/2)                         (1/2)         
       alpha  - 5 alpha - 5      (-alpha (19 alpha + 48))       1      
       -------------------------------------------------------, - alpha
                                     2                          4      
                              4 alpha

You will have to decide which, if any, of these approximations is relevant for your problem.

plot( [eq8], alpha=-2*Pi..2*Pi, view=[DEFAULT,-10..10] );

Hoping this is helpful,

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

It's back, again....

Doug Meade 3088
February 15 2008

Yesterday I, once again, had the old editor. Today, the new one is back.

Very frustrating.

Please reset this, again.

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

axis stretching...

Doug Meade 3088
February 14 2008
0 0
The Maple plot window can be resized dynamically. This changes the aspect ratio of the plot. There a few problems wi th this. First, there is no way to prescribe the dimensions of the plot window. If you stretch the plot and then re-execute the plot command the plot will have the same dimensions as the one it replaces. But, if you delete the plot before re-executing the new plot will have the default dimensions. Another troubling point is that there is no way to precisely set the dimensions of the plot. One consequence of this is that once you change the size of a plot you can never be sure you get it back to its original size. But, maybe this is not what you are wanting to do. Maybe you just want to set the viewing window. This is done with the view= option in any plot command. (Note that display will determine the best viewing window based on the plots that are to be displayed - unless you explicitly put view= in the display command.) Use DEFAULT for any dimensions that you do not want to specify, e.g., view=[DEFAULT, -10..10]. Lastly, I will mention the scaling=constrained optional argument for plot commands. Use this to force the same scale to be used for both the horizontal and vertical dimensions. This is sometimes needed, e.g., to make circles look like circles. I hope one or more of these ideas is appropriate to what you are trying to do. Doug 
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

axis stretching...

Doug Meade 3088
February 14 2008
The Maple plot window can be resized dynamically. This changes the aspect ratio of the plot. There a few problems wi th this. First, there is no way to prescribe the dimensions of the plot window. If you stretch the plot and then re-execute the plot command the plot will have the same dimensions as the one it replaces. But, if you delete the plot before re-executing the new plot will have the default dimensions. Another troubling point is that there is no way to precisely set the dimensions of the plot. One consequence of this is that once you change the size of a plot you can never be sure you get it back to its original size. But, maybe this is not what you are wanting to do. Maybe you just want to set the viewing window. This is done with the view= option in any plot command. (Note that display will determine the best viewing window based on the plots that are to be displayed - unless you explicitly put view= in the display command.) Use DEFAULT for any dimensions that you do not want to specify, e.g., view=[DEFAULT, -10..10]. Lastly, I will mention the scaling=constrained optional argument for plot commands. Use this to force the same scale to be used for both the horizontal and vertical dimensions. This is sometimes needed, e.g., to make circles look like circles. I hope one or more of these ideas is appropriate to what you are trying to do. Doug 
---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

done, but not done...

Doug Meade 3088
February 13 2008

I followed the directions in Will's post, and he replied that he had reset my account so I would not be using the editor. All was fine, or so I thought.

So, why am I seeing the editor as I type this?

To make matters more interesting, when preparing my previous post I was sometimes using the editor and sometimes not. After previewing, most of the time the editor was present but a few times it was not.

It sounds as though something is slightly unstable.

Looking forward to having this resolved.

Doug

After previewing this post, it was completely corrupted. My carefully entered HTML source code was interspersed with copies of the above hyperlink. Strange! And, frustrating!!

---------------------------------------------------------------------

Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/

Agreed!...

Doug Meade 3088
February 08 2008

I'll second acer's comment ...

and re-emphasize that I now type more explicit HTML in my posts than I did with the previous editor.

Doug

---------------------------------------------------------------------
Douglas B. Meade
Math, USC, Columbia, SC 29208  E-mail: mailto:meade@math.sc.edu       
Phone:  (803) 777-6183         URL:    http://www.math.sc.edu/~meade/
First 49 50 51 52 53 54 55 Last Page 51 of 76