Carl Love

Carl Love

28045 Reputation

25 Badges

12 years, 332 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

I just uploaded the worksheet for my numeric solution and put it as a Reply to my numeric Answer below. There is a disturbing difference in the plot of V(t). The symbolic solution is more believable: Since L(t) is nearly linear, ode1 implies that V(t) is nearly constant.

I am going to try adjusting error parameters on the numeric solution.

 



restart:

ode1:=diff(L(t),t)*(V(t))=0.682*10^(-7):

ode2:=137*10^6*diff(V(t),t)*V(t)=(1/(4*3.14))*1.66*10^(-10)*L(t)^2*(3.5*10^3-V(t))^2:

ic1:=V(0)=0.01:

ic2:=L(0)=0:

sys:= {ode1,ode2,ic1,ic2}, {V(t),L(t)}, numeric:

Sol:= dsolve(sys):

plots:-odeplot(Sol, [t,L(t)], t= 0..100);

plots:-odeplot(Sol, [t,V(t)], t= 0..100);

The plot for V(t) is not believable. I expect something nearly constant. But if I narrow the range:

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

Ah ha.

 

Try an implict numeric solver.

Sol:= dsolve(sys, method= mebdfi):

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

Sol:= dsolve(sys, method= lsode):

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

Sol:= dsolve(sys, method= rosenbrock, range= 0..15):

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

 



Download moji_numeric.mw

 



restart:

ode1:=diff(L(t),t)*(V(t))=0.682*10^(-7):

ode2:=137*10^6*diff(V(t),t)*V(t)=(1/(4*3.14))*1.66*10^(-10)*L(t)^2*(3.5*10^3-V(t))^2:

ic1:=V(0)=0.01:

ic2:=L(0)=0:

sys:= {ode1,ode2,ic1,ic2}, {V(t),L(t)}, numeric:

Sol:= dsolve(sys):

plots:-odeplot(Sol, [t,L(t)], t= 0..100);

plots:-odeplot(Sol, [t,V(t)], t= 0..100);

The plot for V(t) is not believable. I expect something nearly constant. But if I narrow the range:

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

Ah ha.

 

Try an implict numeric solver.

Sol:= dsolve(sys, method= mebdfi):

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

Sol:= dsolve(sys, method= lsode):

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

Sol:= dsolve(sys, method= rosenbrock, range= 0..15):

plots:-odeplot(Sol, [t,V(t)], t= 0..15);

 



Download moji_numeric.mw

@N00bstyle 

I'm sorry, I led you astray with SUBS (3). The `if` command has very special evaluation rules that I don't fully understand. Consequently, the subs needs to be done twice. Below, I make the "substituter" a local procedure so that the whole subs command does not need to entered twice.

modalmasspedestriansfail:= proc(_h)
local Ev:= x-> subs(height= _h, x);
     Ev(
           `if`(
                 Ev(qpercentagedifference) >= 5/100
,
              modalmasswithpedestrians,
                 modalmasswithoutpedestrians
           )
     )
end proc:

or

modalmasspedestriansfail:= proc(_h)
local Ev:= x-> subs(height= _h, x);
     `if`(
            Ev(qpercentagedifference) >= 5/100
,
          Ev(modalmasswithpedestrians),
            Ev(modalmasswithoutpedestrians)
     )
end proc:

@N00bstyle 

I'm sorry, I led you astray with SUBS (3). The `if` command has very special evaluation rules that I don't fully understand. Consequently, the subs needs to be done twice. Below, I make the "substituter" a local procedure so that the whole subs command does not need to entered twice.

modalmasspedestriansfail:= proc(_h)
local Ev:= x-> subs(height= _h, x);
     Ev(
           `if`(
                 Ev(qpercentagedifference) >= 5/100
,
              modalmasswithpedestrians,
                 modalmasswithoutpedestrians
           )
     )
end proc:

or

modalmasspedestriansfail:= proc(_h)
local Ev:= x-> subs(height= _h, x);
     `if`(
            Ev(qpercentagedifference) >= 5/100
,
          Ev(modalmasswithpedestrians),
            Ev(modalmasswithoutpedestrians)
     )
end proc:

@N00bstyle 

In method (2), you have the [] inside the parentheses; it needs to go outside, as in this skeleton:

A,B,C:= subs(...= ..., [..., ..., ...] ) []

This says "Make the substitutions in the list [..., ..., ...], THEN turn the list into a sequence (of three items) to match with A,B,C." With the [] inside the parentheses, it says "Turn the list into a sequence, THEN pass that to subs," which subs interprets as nonsense because you are trying to pass it four arguments.

 

In method (3), you used the wrong type of quote marks. It's `if`, not 'if'. On a standard US keyboard, the single back quote is in the upper left corner, under ESC, on the same key as ~.

@N00bstyle 

In method (2), you have the [] inside the parentheses; it needs to go outside, as in this skeleton:

A,B,C:= subs(...= ..., [..., ..., ...] ) []

This says "Make the substitutions in the list [..., ..., ...], THEN turn the list into a sequence (of three items) to match with A,B,C." With the [] inside the parentheses, it says "Turn the list into a sequence, THEN pass that to subs," which subs interprets as nonsense because you are trying to pass it four arguments.

 

In method (3), you used the wrong type of quote marks. It's `if`, not 'if'. On a standard US keyboard, the single back quote is in the upper left corner, under ESC, on the same key as ~.

@erik10 The () makes the procedure immediately above it execute; it's like the () in f(). Without the () it would just be a procedure definition.

@erik10 The () makes the procedure immediately above it execute; it's like the () in f(). Without the () it would just be a procedure definition.

@erik10 I forgot to test my code after a restart. I just had to add two lines at the top, which you'll see. Here's the corrected version.

Simu_para_evalhf.mw

@erik10 I forgot to test my code after a restart. I just had to add two lines at the top, which you'll see. Here's the corrected version.

Simu_para_evalhf.mw

@erik10 Erik wrote:

In a previous comment you suggested to remove my variable broekdel, but I have introduced it again, because it increases the execution time, when the program need to calculate the left hand side in A[number] - j  > L[n] again and again.

Yes, I realized after I posted that removing that variable wouldn't be a good idea if you decided to make c (the number of columns) larger.

Erik wrote:

One last thing I would like is to use the parallelism suggested by Carl, but I am not sure how to combine your earlier code with Acers.

But the last version that I posted does have a parallelization of Acer's code! That was whole point of my latest version.

@erik10 Erik wrote:

In a previous comment you suggested to remove my variable broekdel, but I have introduced it again, because it increases the execution time, when the program need to calculate the left hand side in A[number] - j  > L[n] again and again.

Yes, I realized after I posted that removing that variable wouldn't be a good idea if you decided to make c (the number of columns) larger.

Erik wrote:

One last thing I would like is to use the parallelism suggested by Carl, but I am not sure how to combine your earlier code with Acers.

But the last version that I posted does have a parallelization of Acer's code! That was whole point of my latest version.

Here, I generate the same sequence in Maple. This works nearly instantaneously. But it is still quite impressive that Mathematica figures all this out just from the command that you gave.


restart:

Let  A(n) be the "power tower" sequence, defined recursively by A(1) = 1, A(n) = n^A(n-1). The following procedure computes A(n) mod 10^10, i.e., the last ten decimal digits of A(n).

ph:= numtheory:-phi(10^10):

PowerTower:= proc(n)
option remember;
     n &^ (`if`(n<5, 0, ph)+thisproc(n-1)) mod 10^10
end proc:

PowerTower(1):= 1:

for n to 50 do printf("%2d...%10d\n", n, PowerTower(n)) end do;

 1...         1

 2...         2
 3...         9
 4...    262144
 5...8212890625
 6...1787109376
 7...9058585601
 8...6797099008
 9...8779806721
10...         0
11...         1
12...1445312512
13... 116372481
14...2435774464
15...8212890625
16...1787109376
17... 266721281
18...2305824768
19...5482787841
20...         0
21...         1
22...9316406272
23...7653969921
24...6742273024
25...8212890625
26...1787109376
27...1474856961
28...6950320128
29... 620559361
30...         0
31...         1
32...7187500032
33...6581314561
34...  27920384
35...8212890625
36...1787109376
37...2682992641
38...9565913088
39...4864097281
40...         0
41...         1
42...5058593792
43...7786342401
44... 543340544
45...8212890625
46...1787109376
47...3891128321
48...2131611648
49...8387737601
50...         0

 


Download PowerTower.mw

The above plot, which is surprisingly smooth for an error plot, was made at my default setting of Digits, 15. At the more usual setting Digits = 10, you get the following more-erratic and more-normal error plot. Note that the magnitude of the scale is the same. It is also interesting to see how adjusting the dsolve options abserr, initmesh, and maxmesh affects the plot (see ?dsolve,numeric,BVP ). Note that the default setting of abserr is 1e-6, and the errors that we're actually getting are several orders of magnitude below that. So, that's fairly impressive.

First 642 643 644 645 646 647 648 Last Page 644 of 709