acer

32333 Reputation

29 Badges

19 years, 319 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

A simple way (though likely not the most efficient, if you want to do it thousands of times),

A:=<2,3>;

                                       [2]
                                  A := [ ]
                                       [3]

A . <1,1>;

                                      5

acer

Try changing the call to dsolve(...,numeric) to be,

sol1 := dsolve({DE1, DE2, ICs, a(t) = diff(x(t), t, t)}, numeric, range = 0 .. 3);

Then, your last plot could be produced using,

odeplot(sol1, [t, a(t)], 0 .. N);

That could be a way to get better results, numerically. It makes dsolve(...,numeric) compute a(t) itself. By instead using D(v)(t) then when that gets evalf'd Maple will internally call fdiff() to compute a derivative numerically using a multi-point difference scheme. I would expect dsolve(...,numeric) to be able to compute a(t) more robustly than that, in general if not specifically.

acer

Please check whether these results are plausible and accurate enough, etc.

dsparint2.mw

This revision of code as attached worksheet, including plots out to t=48, takes approximately 61 sec in 64bit Maple 15 on Windows 7 running on an Intel i7 930.

The call to forget(`evalf/int`) is necessary when using operator form for the integrand passed to evalf(Int(..)), if that integrand evaluates according to some hidden quality such as the setting of earlier dsolve/numeric parameters.  

I believe that the call to forget(evalf) may not be necessary and so I removed it. That brought a big performance change (if accompanied by the other edits). I don't know for sure that evalf is not storing any results mistakenly (by ignoring the hidden depencies upon the parameters of the earlier dsolve/numeric output procs). I didn't have the patience to wait the long time to get final results using your posted code.

I changed it to use evalf(Int(...)) instead of evalf(int(...)). But to do that, with the forced method=_d01ajc, I had to workaround what appears to be a bug in evalf/Int. It looks like an internal evalf/Int routine might be unprepared to catch & handle a particular error coming out of evalhf when trying to run the earlier dsolve/numeric stuff during evaluation of the integrand. And so it doesn't seem able to gracefully fall back to regular evalf mode. I (believe) that I worked around this by the unusual maneouver of turning the integrand into a procedure which is deliberately non-evalhf'able (viz. it creates an empty list) and which would throw a different but recognized evalhf error. Blech.

And I made a few adjustments to try and tone down some of the accuracy tolerances.

nb. I only changed Po2m to something usable, since Pco2m doesn't seem needed by the final dsolve.

The final plot of O2r looks a bit different in Maple 16 than it does in Maple 15, especially in the key t=15..25 range. The other four plotted functions of t looked pretty much the same in both versions. But in Maple 16 O2r doesn't really fall below 0.00004 whereas in Maple 15 it gets down to -0.0004 or so. I didn't try another solver for dsolve/numeric. But this discrepency seems to persist even if I adjust the final dsolve relerr to be 1e-9 and the evalf/Int epsilon to be 1e-11.

acer

w:=-t*sin(x*Pi): # or what have you...

A:=plots:-animate(w, x=0..1, t=0..0.1, frames=50, numpoints=100,
                  color=red, thickness=3, view=[0..1,-1..1]):

plots:-display(A, view=[0..1,-1..1], insequence=true);

plots:-display(plottools:-transform((x,y)->[y,x])(A),
               view=[-1..1,0..1], insequence=true);

acer

See the help-page evalf,Int which is also cross-referenced from the 5th bullet point in the Description section of the help-page for int.

For some other  examples you could look again at a response to your question from last week.

acer

In Maple 16.00 running on Windows 7 with 6GB RAM available, I get results as follows from kernelopts(bytesalloc),  for a single set of attempts.

32bit Maple: n=100  171MB

64bit Maple: n=100  264MB

acer

You could use plots:-odeplot,

gamma1:=1:
gamma2:=1:

l:=t->sqrt((x1(t)-x2(t))**2+(y1(t)-y2(t))**2):

DE:=[diff(x1(t),t)=(-1/(2*Pi))*gamma2*(y1(t)-y2(t))/l(t),
     diff(y1(t),t)=(-1/(2*Pi))*gamma1*(x1(t)-x2(t))/l(t),
     diff(x2(t),t)=(-1/(2*Pi))*gamma2*(y2(t)-y1(t))/l(t),
     diff(y2(t),t)=(-1/(2*Pi))*gamma1*(x2(t)-x1(t))/l(t)]:

plots:-odeplot(dsolve({DE[],x1(0)=1,x2(0)=-1,y1(0)=1,y2(0)=0},numeric),
               [[x1(t),y1(t)],[x2(t),y2(t)]],t=0..50,frames=25,
               color=[gold,gold],thickness=2);

acer

Don't break the long name that looks like `#mrow...`

with(Statistics):
BoxPlot([seq(Sample(Normal(ln(i), 3), 100), i = 1..20)],
title = `#mrow(mi("Jan 2011",mathcolor="#ff00ff"),mi(" --- "),mi("Dec 2011",mathcolor="#005555"))`,
color="ff00ff".."005555",deciles = false);

Statistics:-BoxPlot([seq(Statistics:-Sample(Normal(ln(i), 3), 100), i = 1..20)], title = `#mrow(mi(

Another example, building up that long TypeMK name (which you could also do with a procedure),

clr1:="#999900":
clr2:="#009999":

Statistics:-BoxPlot(
    [seq(Statistics:-Sample(Normal(ln(i),3),100),i=1..10)],
    title = cat(`#mrow(`,
            `mi("Jan 2011",mathcolor=\"`,clr1,`\")`,
            `,`,
            `mi(" --- ")`,
            `,`,
            `mi("Dec 2011",mathcolor=\"`,clr2,`\")`,
             `)`),
    color = clr1..clr2, deciles = false);

Statistics:-BoxPlot(     [seq(Statistics:-Sample(Normal(ln(i),3),100),i=1..10)],     title = cat(`#mrow(`,             `mi("Jan 2011",mathcolor=\"`,"#999900",`\")`,             `,`,             `mi(" --- ")`,             `,`,             `mi("Dec 2011",mathcolor=\"`,"#009999",`\")`,              `)`),     color = "#999900".."#009999", deciles = false)

acer

The explanation is confusing. Did you mean to say instead that the entries of N1, N2... are expressions in terms of `eta`?

Are you trying to map the `int` or `Int` commands over N1,N2... with eta=eta0..eta2 being the extra argument? If so that the syntax is not right. It should be something like map(int,N1,eta=eta0..eta2) or map(Int,N1,eta=eta0..eta2), in that case.

Do you have somethin like this, but for 3x3 Matrices,

N1:=Matrix([[eta,eta^2],[sin(eta),cos(eta)*eta]]);

                             [                 2    ]
                             [  eta         eta     ]
                       N1 := [                      ]
                             [sin(eta)  cos(eta) eta]

L:=[eta0=8.3,eta2=666];

                        L := [eta0 = 8.3, eta2 = 666]

# First way: inert integral with parameters substituted by
# values, then numerical integration performed.
# Be careful with syntax. And note Int not int.
E1:=map(Int,N1,eta=eta0..eta2):
eval(E1,L):
evalf(%);

                     [              5                7]
                     [2.217435550 10   9.846924140 10 ]
                     [                                ]
                     [   -1.431221219     -17.80614111]

# Alternatively, exact integration attempted, followed by
# floating-point evaluation of that result.
E1:=map(int,N1,eta=eta0..eta2):
eval(E1,L):
evalf(%);

                     [              5                7]
                     [2.217435550 10   9.846924140 10 ]
                     [                                ]
                     [   -1.431221219     -17.80614111]

In both ways, you don't have to assign to `eta0` or `eta2`, or change the names used. You just have to build the various lists L.

acer

You have an E(T) which should likely be E(t). You need to use the infix `union` command to merge the sets `dsys` and `isc`.

I changed e^7 to exp(7) which how one exponentiates the base of the natural logarithm in Maple. Or did you intend to use `e` as a variable parameter?

restart;

dsys := {diff(E(t), t) = -kf*E(t)*S(t)+kr*ES(t)+kcat*ES(t),

                diff(ES(t), t) = kf*E(t)*S(t)-kr*ES(t)-kcat*ES(t),

                diff(P(t), t) = kcat*ES(t),

                diff(S(t), t) = -kf*E(t)*S(t)+kr*ES(t)};
  
kf := 10^6; kr := 10^(-4); kcat := .1;
                        
isc := {E(0) = 2/exp(7), ES(0) = 0, P(0) = 0, S(0) = 5/exp(7)};
        
dsol := dsolve(dsys union isc, numeric, method = gear[polyextr],
                initstep = 0.15e-1, minstep = 10^(-11),
                abserr = 0.1e-4, relerr = 0.1e-5);

acer

It's a square.

factor(dell);

                                         2
                      /          (1/2)  \ 
                      \-a + 9 + 6      b/ 

acer

You have X with 8 entries and Y with 9 entries.

Does it help if you reduce the number of entries of Y by one, so as to match that of X?

acer

M := Matrix([[1,2],[3,4],[5,6]]);


                                      [1  2]
                                      [    ]
                                 M := [3  4]
                                      [    ]
                                      [5  6]

M . Vector([1/2,1/2]);

                                    [ 3]
                                    [ -]
                                    [ 2]
                                    [  ]
                                    [ 7]
                                    [ -]
                                    [ 2]
                                    [  ]
                                    [11]
                                    [--]
                                    [2 ]

M . Vector([0.5,0.5]);

                             [1.50000000000000]
                             [                ]
                             [3.50000000000000]
                             [                ]
                             [5.50000000000000]

acer

Here is a version that runs in Maple 12.

The difficulty is the quickly oscillating semi-infinite integral involving BesselJ(1,..) which is decaying slowly. I split the range from 0..infinity to 0..1 and 1..infinity, and then changed variables on the second of those. That gets a sum of two oscillating integrals on 0..1, and these can be handled by method=_d01akc.

It's not particularly fast to compute. The plot takes about a minute on a fast Intel i7, because each evaluation of func is not very fast. I made func a procedure (both to try and get some more speed when doing the outer integral after the plot, but also to try and control the working precision forcibly).

The inner integral involving the Bessel J1 needs to be computed accurately enough that its results allow the final outer integral to converge.

 

restart:

func:=proc(r) option remember;
  .2359951354e-4/Pi*exp(-800000000.0*(r-.375e-4)^2) *
  (evalf(Int(k->-k^2*BesselJ(1.,k*r)*exp(-.1757812500e-9*k^2)/(.66e-2*k^2+2423.07),
             0. .. 1, digits=14, epsilon=1e-11, method=_d01akc)) +
   evalf(Int(K->-BesselJ(1., r/K)*exp(-0.1757812500e-9/K^2)/(K^4*(0.66e-2/K^2+2423.07)),
             0. .. 1., digits=14, epsilon=1e-11, method=_d01akc)));
end proc:

 

st:=time():
plot(func, -0.200e-3 .. 0.200e-3); #, adaptive=false, numpoints=20);
time()-st;

 

78.547

(1)

st:=time():
evalf(Int(func,-0.100e-3 .. 0.100e-3,digits=14,epsilon=1e-7,method=_d01ajc));
time()-st;

-0.11785650145932e-2

 

54.881

(2)

 

 

Download inty.mw

acer

First 265 266 267 268 269 270 271 Last Page 267 of 336