acer

32385 Reputation

29 Badges

19 years, 334 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Here is another way to handle this example,

restart;

kernelopts(version);

   Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874

fprime := x-> x^6-(3/2)*x^5+2*x^4+(5/2)*x^3-7*x^2+2:
f := unapply(simplify(int(fprime(x), x)), x):
g := unapply(expand(f(x^2+2*x)), x):

sols := [solve({x>-infinity, x<infinity,
                 diff(g(x),x)=0, diff(g(x),x,x)<>0},x)]:

evalf(sols);

      [{x = -1.}, {x = 0.2834077704}, {x = -2.283407770}, 
       {x = 0.4977473282}, {x = -2.497747328}, {x = -0.3051050611}, 
       {x = -1.694894939}]

sort(map[2](eval,x,evalf(sols)));

    [-2.497747328, -2.283407770, -1.694894939, -1., -0.3051050611, 
      0.2834077704, 0.4977473282]

I was able to compile the C file and link the shared library, in a Linux terminal.

Using the commandline interface for Maple 2019.0 in that Linux terminal I was get it to run from within Maple after a call to define_external.

I added a few lines in the C source to printf the entries, before the computational portion. And I threw in a printf("\n") before the return(0).

Here is what my terminal session looked like:

restart:                                                                                     

loc:=cat(kernelopts(homedir),"/gauss.so"):                                                   

mygauss:=define_external('main',RETURN::integer[4],LIB=loc):                                 

mygauss();                                                                                   

Enter the order of matrix: 
Enter the elements of augmented matrix row-wise:

A[1][1] : A[1][2] : 
A[1][1] = 2.000000

A[1][2] = 3.000000

The solution is: 

x1=1.500000	
                          0

As Carl has mentioned, using scanf in C to input the entries is unusual. It's awkward, and possibly not acting just as intended (prematurely accepting and receiving input values, say). It behaved a little weird for me -- I'm not sure that I'd easily be able to enter a 2x2 example with float entries. You'd be far better off specifying Matrix arguments in the define_external, both to pass in the data, and to act in-place to hold the result.

But the first issue is whether the define_external will succeed for you.

Are you 100% sure that the LIB location that you're using actually points to the correct location of the shared library? Is your directory actually called "/home/user/Desktop"? Where is the shared library file? Is it somewhere like, say, "/home/radaar/gauss.so" ?

Is this acceptable?

I wrote it to try and handle only the kind of indexing in your example, eg, C[1][i,m] etc. No doubt something recursive could be concocted to handle arbitrary depth of indexing and arbitrary degree of indexing (per depth).

restart;

ee := sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M)-A__0*H__1*(sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M))-A__1*F*(sum(C[2][i, k]*psi[theta][k, j], k = 1 .. N))+A__1*K__4*(sum(C[1][j, k]*psi[theta][k, i], k = 1 .. N))/r__i+A__0*K__4*(sum(C[1][j, k]*v[0][k, i], k = 1 .. N))/(2*r__i)

sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M)-A__0*H__1*(sum(sum(C[1][i, m]*C[1][j, k]*v[0][k, m], k = 1 .. N), m = 1 .. M))-A__1*F*(sum(C[2][i, k]*psi[theta][k, j], k = 1 .. N))+A__1*K__4*(sum(C[1][j, k]*psi[theta][k, i], k = 1 .. N))/r__i+(1/2)*A__0*K__4*(sum(C[1][j, k]*v[0][k, i], k = 1 .. N))/r__i

new:=subsindets(ee,
                And(indexed,satisfies(u->op(0,u)::indexed)),
                u->nprintf(`%a__%a__%s`,
                           op([0,0],u),op([0,..],u),
                           cat(op(1,u),
                               seq(cat(`,`,p),p=[op(2..,u)]))));

sum(sum(`C__1__i,m`*`C__1__j,k`*`v__0__k,m`, k = 1 .. N), m = 1 .. M)-A__0*H__1*(sum(sum(`C__1__i,m`*`C__1__j,k`*`v__0__k,m`, k = 1 .. N), m = 1 .. M))-A__1*F*(sum(`C__2__i,k`*`psi__theta__k,j`, k = 1 .. N))+A__1*K__4*(sum(`C__1__j,k`*`psi__theta__k,i`, k = 1 .. N))/r__i+(1/2)*A__0*K__4*(sum(`C__1__j,k`*`v__0__k,i`, k = 1 .. N))/r__i

latex(eval(new,1));

\sum _{m=1}^{M} \left( \sum _{k=1}^{N}C_{\mbox {{\tt 1}}_{\mbox {{\tt

i,m}}}}\,C_{\mbox {{\tt 1}}_{\mbox {{\tt j,k}}}}\,v_{\mbox {{\tt 0}}_{
\mbox {{\tt k,m}}}} \right) -A_{0}\,H_{1}\,\sum _{m=1}^{M} \left(
\sum _{k=1}^{N}C_{\mbox {{\tt 1}}_{\mbox {{\tt i,m}}}}\,C_{\mbox {{
\tt 1}}_{\mbox {{\tt j,k}}}}\,v_{\mbox {{\tt 0}}_{\mbox {{\tt k,m}}}}
 \right) -A_{1}\,F\sum _{k=1}^{N}C_{\mbox {{\tt 2}}_{\mbox {{\tt i,k}}
}}\,\psi_{\theta_{\mbox {{\tt k,j}}}}+{\frac {A_{1}\,K_{4}\,\sum _{k=1
}^{N}C_{\mbox {{\tt 1}}_{\mbox {{\tt j,k}}}}\,\psi_{\theta_{\mbox {{
\tt k,i}}}}}{r_{i}}}+1/2\,{\frac {A_{0}\,K_{4}\,\sum _{k=1}^{N}C_{
\mbox {{\tt 1}}_{\mbox {{\tt j,k}}}}\,v_{\mbox {{\tt 0}}_{\mbox {{\tt
k,i}}}}}{r_{i}}}

 

Download end_ac.mw

seq(x[i], i = 1 .. 4);
 
                x[1], x[2], x[3], x[4]

seq(cat(x,i), i = 1 .. 4);

                     x1, x2, x3, x4

Here are the plots of Re(c1) and Im(c1).

I used two different formulas for A and S1 (since above you have used ones which do not correspond to the correct translation from Mathematica, that you asked about yesterday).

It looks as if the revised formulas produce results with zero imaginary component (and so might be what you were expecting?).

restart:

kernelopts(version);

`Maple 13.00, X86 64 LINUX, Apr 13 2009, Build ID 397624`

(1)

Digits:=15:

d1:=0.2:L1:=0.2:L2:=0.2:B1:=0.7:B:=1:beta:=0.01:

d2:=0.6:m:=0.1:k:=0.1:

h:=z->piecewise( z<=d1,    1,

                z<=d1+L1,   1-(gamma1/(2))*(1 + cos(2*(Pi/L1)*(z-d1-L1/2))),

                z<=B1-L2/2,  1 ,          

                z<=B1,  1-(gamma2/(2))*(1 + cos(2*(Pi/L2)*(z - B1))),

                z<=B1+L2/2,  1-(gamma2/(2))*(1 + cos(2*(Pi/L2)*(z - B1))),

               z<=B,    1):

A:=(-m^2/4)-(1/(4*k)):

S1:=(h(z)^2)/(4*A)-ln(A*h(z)^2+1)*(1+h(z)^2)/(4*A):

b1:=1/S1:

c1:=Int(b1,z=0..1,method=_Dexp,digits=15,epsilon=1e-4):

plot([seq(eval(Re(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Re(c1)'],
     gridlines=false, axes=boxed);

 

plot([seq(eval(Im(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Im(c1)'],
     gridlines=false, axes=boxed);

 

A:=(-m^2/4)-(1/4*k):

S1:=(h(z)^2)/4*A-ln(A*h(z)^2+1)*(1+h(z)^2)/4*A:

b1:=1/S1:

c1:=Int(b1,z=0..1,method=_Dexp,digits=15,epsilon=1e-4):

plot([seq(eval(Re(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Re(c1)'],
     gridlines=false, axes=boxed);

 

plot([seq(eval(Im(c1),gamma2=j),j in[0,0.02,0.06])],gamma1=0.02..0.1,
     legend = [gamma2 = 0.0,gamma2 = 0.02,gamma2 = 0.04],
     linestyle = [solid,dash,dot],
     color = [black,black,black],
     labels=[gamma1,'Im(c1)'],
     gridlines=false, axes=boxed);

 

 

Download waseem_plot2.mw

These two approaches work in my Maple 2018.2,

restart;

kernelopts(version);

`Maple 2018.2, X86 64 LINUX, Nov 16 2018, Build ID 1362973`

vals:=[k=1.3806E-23, T=300, h__bar=1.05456E-34,
       L__z=6E-9, m__hhw=3.862216E-31,
       E__hh0=3.012136E-21, E__hh1=1.185628E-20]:

p := k*T*m__hhw*(ln(1+exp(E__fv-E__hh0)/(k*T))
                 +ln(1+exp(E__fv-E__hh1)/(k*T)))/(Pi*h__bar^2*L__z)
     = 0.3e25; #convert(0.3e25, rational,exact);

k*T*m__hhw*(ln(1+exp(E__fv-E__hh0)/(k*T))+ln(1+exp(E__fv-E__hh1)/(k*T)))/(Pi*h__bar^2*L__z) = 0.3e25

[solve({combine(p), E__fv<0},{E__fv})] assuming E__fv::real, op(vals);
eval(%,vals);

[{E__fv = E__hh0+ln(.5000000000*(-1.*exp(E__hh0-1.*E__hh1)*T*k-1.*k*T+((exp(E__hh0-1.*E__hh1))^2*T^2*k^2+4.*exp(E__hh0-1.*E__hh1)*exp(0.9424777962e25*h__bar^2*L__z/(k*T*m__hhw))*T^2*k^2-2.*exp(E__hh0-1.*E__hh1)*T^2*k^2+k^2*T^2)^(1/2))/exp(E__hh0-1.*E__hh1))}]

[{E__fv = -48.46001884}]

subsindets(combine(p),specfunc(ln),u->(ln(expand(op(u))))) assuming E__fv::real, op(vals);
eval(%,vals);
solve({%, E__fv<0}, {E__fv});

k*T*m__hhw*ln(1+exp(E__fv)/(exp(E__hh1)*k*T)+exp(E__fv)/(exp(E__hh0)*k*T)+(exp(E__fv))^2/(exp(E__hh0)*k^2*T^2*exp(E__hh1)))/(Pi*h__bar^2*L__z) = 0.3e25

0.7631009085e25*ln(1+0.4828818388e21*exp(E__fv)+0.5829371757e41*(exp(E__fv))^2) = 0.3e25

{E__fv = -48.46001884}

 

Download solve_M2018.2.mw

Note that "1/4 k" in Mathematica translates to "1/4*k" in Maple, but you have it wrongly as "1/(4*k)" in Maple. You made a similar mistake in at least two other places.

restart:
epsilon := 0.2:
h := z->1+epsilon*sin(2*Pi*z):
m := 10: k := 0.8:
A := (-m^2/4)-(1/4*k):
S1 := z->(h(z)^2)/4*A-ln(A*h(z)^2+1)*(1+h(z)^2)/4*A:

S1(0.9);

         27.86472748 + 35.20420042 I

Perhaps it's the case that your s or t (or both) can be taken to be real-valued.

restart;

kernelopts(version);

`Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874`

with(LinearAlgebra):

N26 := Vector[row]([ 2*s^2*cos(t), 2*s^2*sin(t), s ]);

Vector[row](3, {(1) = 2*s^2*cos(t), (2) = 2*s^2*sin(t), (3) = s})

combine(Norm(N26, 2)) assuming t::real;

(4*abs(s)^4+abs(s)^2)^(1/2)

simplify(Norm(N26, 2)) assuming t::real;

abs(s)*(4*abs(s)^2+1)^(1/2)

combine(Norm(N26, 2)) assuming real;

(4*s^4+s^2)^(1/2)

simplify(Norm(N26, 2)) assuming real;

(4*s^2+1)^(1/2)*abs(s)

 

Download simp_real.mw

See also simplify(Norm(N26, 2, conjugate=false)) 

restart;

kernelopts(version);

  Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874

f:=ln(s + 2)^2 + 2*polylog(2, -1 - s) + 2*polylog(2, (1 + s)/(s + 2)):

simplify(convert(f,dilog));

                      0
Related examples have been reported previously.
restart;

kernelopts(version);

   Maple 2019.1, X86 64 LINUX, May 21 2019, Build ID 1399874

u := < -(1/4)*a, -(1/12)*sqrt(3)*a, -x >:
v := < -(1/2)*a, (1/6)*sqrt(3)*a, (1/2)*x >:
Auv := Student:-MultivariateCalculus:-Angle(u, v):

solve({Auv = arccos(2*sqrt(26)*(1/13)), x>0, a>0}, x) assuming a>0;

                       /    1    (1/2)  \ 
                      { x = - 114      a }
                       \    6           / 
The inequality a>0 passed to `solve` may not actually be used of necessary. And the assumption on `a` may just serve to simplify the piecewise parametric result from `solve`. It's sad that `solve'` functionality is so disjointed (and no longer properly documented).

Are you supposing that tau and sigma1 are real?

restart;

SS1 := -4*sqrt(sigma1^2+4*tau^2)*tau^2/((sigma1+sqrt(sigma1^2
       +4*tau^2))^2*(1+4*abs(tau/(sigma1+sqrt(sigma1^2+4*tau^2)))^2));

-4*(sigma1^2+4*tau^2)^(1/2)*tau^2/((sigma1+(sigma1^2+4*tau^2)^(1/2))^2*(1+4*abs(tau/(sigma1+(sigma1^2+4*tau^2)^(1/2)))^2))

SS2 := 4*sqrt(sigma1^2+4*tau^2)*tau^2/((-sigma1+sqrt(sigma1^2
       +4*tau^2))^2*(1+4*abs(tau/(-sigma1+sqrt(sigma1^2+4*tau^2)))^2));

4*(sigma1^2+4*tau^2)^(1/2)*tau^2/((-sigma1+(sigma1^2+4*tau^2)^(1/2))^2*(1+4*abs(tau/(-sigma1+(sigma1^2+4*tau^2)^(1/2)))^2))

K1:=simplify(SS1) assuming real;

-2*tau^2*(sigma1^2+4*tau^2)^(1/2)/(sigma1*(sigma1^2+4*tau^2)^(1/2)+sigma1^2+4*tau^2)

K2:=simplify(SS2) assuming real;

2*tau^2*(sigma1^2+4*tau^2)^(1/2)/(sigma1^2-sigma1*(sigma1^2+4*tau^2)^(1/2)+4*tau^2)

Q:=[solve({S1=K1,S2=K2},{tau,sigma1},explicit)];

[{sigma1 = S1+S2, tau = (-S1*S2)^(1/2)}, {sigma1 = S1+S2, tau = -(-S1*S2)^(1/2)}]

diff(sigma1(S1,S2),S1)=diff(eval(sigma1,Q[1]),S1);
diff(sigma1(S1,S2),S2)=diff(eval(sigma1,Q[1]),S2);

diff(sigma1(S1, S2), S1) = 1

diff(sigma1(S1, S2), S2) = 1

diff(tau(S1,S2),S1)=diff(eval(tau,Q[1]),S1);
diff(tau(S1,S2),S2)=diff(eval(tau,Q[1]),S2);

diff(tau(S1, S2), S1) = -(1/2)*S2/(-S1*S2)^(1/2)

diff(tau(S1, S2), S2) = -(1/2)*S1/(-S1*S2)^(1/2)

diff(tau(S1,S2),S1)=diff(eval(tau,Q[2]),S1);
diff(tau(S1,S2),S2)=diff(eval(tau,Q[2]),S2);

diff(tau(S1, S2), S1) = (1/2)*S2/(-S1*S2)^(1/2)

diff(tau(S1, S2), S2) = (1/2)*S1/(-S1*S2)^(1/2)

factor(simplify(SS1+SS2)) assuming real;

sigma1

factor(simplify(SS1*SS2)) assuming real;

-tau^2

 

Download solve_thing.mw

restart;

L1 := [[1,2,3],[7,8,9],[13,12,11]]:

map[2](op,3,L1);
                           [3, 9, 11]

map(`?[]`,L1,[3]);
                           [3, 9, 11]

L2 := [1, [2, 3], [4, [5, 6], 7], [8, 3], 9]:

map[2](op,1,L2);
                        [1, 2, 4, 8, 9]

map(u->`if`(u::list,u[1],u),L2);
                        [1, 2, 4, 8, 9]

L3 := [[1,2,3],[7,8,2],[13,12,1]]:

sort(L3, (a,b) -> a[3]<b[3]);
              [[13, 12, 1], [7, 8, 2], [1, 2, 3]]

If your list contains a mix of list and scalars (like L2 above) then don't use that op approach if the entries are not all lists or integers. And you may also want to allow for empty sublists. For example,

L4 := [1,[2,3],11.02,[],[4,[5,6],7],[8,3],9.0003,[]]:

map(u->`if`(u::list,`if`(nops(u)>0,u[1],NULL),u), L4);
          [1, 2, 11.02, 4, 8, 9.0003]

restart

with(plots):

with(CurveFitting):

Digits := 50:

Direct := diff(C(x, t), t) = -.5*(diff(C(x, t), x))+10*(diff(C(x, t), x, x)):

g := proc (t) options operator, arrow; piecewise(0 < t and t <= 2000, 0, 2000 < t and t <= 4000, 0.5e-2*t-10, 4000 < t and t <= 6000, (-1)*0.25e-2*t+20, 6000 < t and t <= 9000, 0.5e-2*t-25, 9000 < t and t <= 16000, (-1)*0.286e-2*t+45.76, 16000 < t and t <= 20000, 0) end proc:

plot(g(t), t = 0 .. 20000):

``

IBCD := C(x, 0) = 0, C(0, t) = g(t), (D[1](C))(10000, t) = 0:

pdsD := pdsolve(Direct, [IBCD], time = t, range = 0 .. 10000, timestep = 60, numeric, spacestep = 25)

module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; end module

NULL

NULL

pts := [Vector[row](100, {(1) = 0, (2) = 101.010101, (3) = 202.020202, (4) = 303.030303, (5) = 404.040404, (6) = 505.0505051, (7) = 606.0606061, (8) = 707.0707071, (9) = 808.0808081, (10) = 909.0909091, (11) = 1010.10101, (12) = 1111.111111, (13) = 1212.121212, (14) = 1313.131313, (15) = 1414.141414, (16) = 1515.151515, (17) = 1616.161616, (18) = 1717.171717, (19) = 1818.181818, (20) = 1919.191919, (21) = 2020.20202, (22) = 2121.212121, (23) = 2222.222222, (24) = 2323.232323, (25) = 2424.242424, (26) = 2525.252525, (27) = 2626.262626, (28) = 2727.272727, (29) = 2828.282828, (30) = 2929.292929, (31) = 3030.30303, (32) = 3131.313131, (33) = 3232.323232, (34) = 3333.333333, (35) = 3434.343434, (36) = 3535.353535, (37) = 3636.363636, (38) = 3737.373737, (39) = 3838.383838, (40) = 3939.393939, (41) = 4040.40404, (42) = 4141.414141, (43) = 4242.424242, (44) = 4343.434343, (45) = 4444.444444, (46) = 4545.454545, (47) = 4646.464646, (48) = 4747.474747, (49) = 4848.484848, (50) = 4949.494949, (51) = 5050.505051, (52) = 5151.515152, (53) = 5252.525253, (54) = 5353.535354, (55) = 5454.545455, (56) = 5555.555556, (57) = 5656.565657, (58) = 5757.575758, (59) = 5858.585859, (60) = 5959.59596, (61) = 6060.606061, (62) = 6161.616162, (63) = 6262.626263, (64) = 6363.636364, (65) = 6464.646465, (66) = 6565.656566, (67) = 6666.666667, (68) = 6767.676768, (69) = 6868.686869, (70) = 6969.69697, (71) = 7070.707071, (72) = 7171.717172, (73) = 7272.727273, (74) = 7373.737374, (75) = 7474.747475, (76) = 7575.757576, (77) = 7676.767677, (78) = 7777.777778, (79) = 7878.787879, (80) = 7979.79798, (81) = 8080.808081, (82) = 8181.818182, (83) = 8282.828283, (84) = 8383.838384, (85) = 8484.848485, (86) = 8585.858586, (87) = 8686.868687, (88) = 8787.878788, (89) = 8888.888889, (90) = 8989.89899, (91) = 9090.909091, (92) = 9191.919192, (93) = 9292.929293, (94) = 9393.939394, (95) = 9494.949495, (96) = 9595.959596, (97) = 9696.969697, (98) = 9797.979798, (99) = 9898.989899, (100) = 10000}), Vector[row](100, {(1) = 0., (2) = 0., (3) = 0., (4) = 0., (5) = 0., (6) = 0., (7) = 0.2e-4, (8) = 0.10e-3, (9) = 0.56e-3, (10) = 0.248e-2, (11) = 0.919e-2, (12) = 0.2857e-1, (13) = 0.7551e-1, (14) = .17180, (15) = .34128, (16) = .60125, (17) = .95511, (18) = 1.39139, (19) = 1.88934, (20) = 2.42694, (21) = 2.98667, (22) = 3.55718, (23) = 4.13232, (24) = 4.70924, (25) = 5.28676, (26) = 5.86447, (27) = 6.44223, (28) = 7.02000, (29) = 7.59778, (30) = 8.17555, (31) = 8.75333, (32) = 9.33110, (33) = 9.90886, (34) = 10.48658, (35) = 11.06420, (36) = 11.64158, (37) = 12.21841, (38) = 12.79408, (39) = 13.36737, (40) = 13.93611, (41) = 14.49665, (42) = 15.04318, (43) = 15.56715, (44) = 16.05673, (45) = 16.49676, (46) = 16.86924, (47) = 17.15449, (48) = 17.33308, (49) = 17.38802, (50) = 17.30713, (51) = 17.08478, (52) = 16.72287, (53) = 16.23085, (54) = 15.62469, (55) = 14.92534, (56) = 14.15692, (57) = 13.34508, (58) = 12.51563, (59) = 11.69367, (60) = 10.90282, (61) = 10.16470, (62) = 9.49811, (63) = 8.91819, (64) = 8.43537, (65) = 8.05441, (66) = 7.77367, (67) = 7.58484, (68) = 7.47331, (69) = 7.41912, (70) = 7.39856, (71) = 7.38612, (72) = 7.35668, (73) = 7.28756, (74) = 7.16043, (75) = 6.96260, (76) = 6.68785, (77) = 6.33655, (78) = 5.91524, (79) = 5.43562, (80) = 4.91320, (81) = 4.36573, (82) = 3.81158, (83) = 3.26826, (84) = 2.75118, (85) = 2.27275, (86) = 1.84190, (87) = 1.46393, (88) = 1.14072, (89) = .87118, (90) = .65192, (91) = .47788, (92) = .34306, (93) = .24113, (94) = .16592, (95) = .11174, (96) = 0.7364e-1, (97) = 0.4749e-1, (98) = 0.2996e-1, (99) = 0.1850e-1, (100) = 0.1244e-1})]

A := unapply(Spline(pts[1], pts[2], x, degree = 2), x):

plot(A(x), x = 0 .. 10000):

Inverse := diff(C(x, t), t) = .5*(diff(C(x, t), x))-50000*(diff(C(x, t), x, x, x, x))-10*(diff(C(x, t), x, x)):

NULL

NULL

IBCI := C(x, 0) = A(x), (D[1](C))(0, t) = 0, (D[1](C))(10000, t) = 0, (D[1, 1](C))(0, t) = 0, (D[1, 1](C))(10000, t) = 0:

pdsI := pdsolve(Inverse, [IBCI], time = t, range = 0 .. 10000, timestep = .5, numeric, spacestep = 25)

module () local INFO; export plot, plot3d, animate, value, settings; option `Copyright (c) 2001 by Waterloo Maple Inc. All rights reserved.`; end module

A1 := pdsI:-animate(t = 20000, frames = 20, title = "time = %f"):

A1;

A2 := pdsD:-animate(t = 20000, frames = 20, title = "time = %f"):

A2;

h := display(A1, A2):

h;

 

Download both_of_them_ac.mw

@ajfriedlan You could have a look at the Grid package (and perhaps Grid:-Map in particular, to map your integrating process across a list/Array of parameter values). The big difference being that the Grid package parallelizes by launching new and separate kernel engines to do the work -- thus potentially alleviating the need for thread-safety.

As for ArrayTools:-Concatenate, I am not sure about that particular procedure but I notice that it calls ArrayTools:-Copy. Now, ArrayTools:-Copy does an interesting thing: the first time it gets called (per restart) it replaces itself with a call to a compiled external function. But that replacement operation is not thread-safe, and if one thread tries to call the command while another is in the process of making its first call then bad things could happen. So one simplistic way to try and deal with this is to execute a (typical, for your code) preliminary, dummy call to ArrayTools:-Concatenate and/or ArrayTools:-Copy before you try can run your code under the Threads package. But if you manage to replace all your Threads attempts with Grid attempts then this all may be moot to you.

 

 

Here's an alternative, which keeps construction of A1 and A2 separate (if that's what you prefer).

The key thing is that the legend becomes a so-called local property of the curve plotting structure (technically, a substructure of the CURVES plotting structure).

That's the case in each frame, with the code below. The plot command knows how to do that. But when you originally called plots:-display([A1,A2],legend=[...],...) then the plots:-display command did not know how to dig into the A1 and A2 animations and find which curves to associate with the supplied legend items.

Sure, that is a weakness of plots:-display (but there are much harder situations in similar examples, and it would be inconsistent to only handle a smattering of cases). 

c := x -> piecewise(0 <= x and x <= 450000, 0.37*x, 450000 < x,
                    0.37*x + 0.06*(x - 450000)):
g := x -> piecewise(0 <= x and x <= 558000, 0.37*x, 558000 < x,
                    0.37*x + 0.12*(x - 558000)):

A1 := plots:-animate(plot, [g(x), x = 0 .. skat, color = blue, legend="Line 1"],
                     frames = 20, skat = 0 .. 1000000):
A2 := plots:-animate(plot, [c(x), x = 0 .. skat, color = red, legend="Line 2"],
                     frames = 20, skat = 0 .. 1000000):

plots:-display([A1, A2], size = [500, 350], gridlines = true,
               legendstyle=[font=[Lucida, roman, 14], location=bottom]);

First 150 151 152 153 154 155 156 Last Page 152 of 336