Advanced Search

Rouben Rostamian

7711 Reputation

22 Badges

19 years, 3 days
Contact Rouben Rostamian

MaplePrimes Activity


These are answers submitted by Rouben Rostamian

Workaround...

Rouben Rostamian 7711
June 16 2018
2 3

That looks like a bug to me.  Maple 2018 does not crash, but fails with the error message "Error, (in Explore) while processing result".

Here is a workaround that works on Maple 2018.  See if it works for you.  Change
caption = typeset(x^2 + y^2)
to
caption = typeset("%1+%2", x^2, y^2)
 

 

typeset...

Rouben Rostamian 7711
June 15 2018
1 0
> 

restart:
> 

omega0:= 10/3.:

Choose one of  legend_plain  or  legend_fancy:

> 

legend_plain := [typeset(omega[0] = omega0)];

[typeset(omega[0] = 3.333333333)]

> 

legend_fancy := [typeset(omega[0] = convert(sprintf("%.1f", omega0), symbol))];

[typeset(omega[0] = `3.3`)]

> 

plot(omega^2/omega0, omega = 1..1.5, legend = legend_fancy);

 

Download mw.mw

Here is how...

Rouben Rostamian 7711
June 11 2018
2 4

I can't tell by looking in your worksheet what you want your animation to show.

But it appears that you know how to produce the individual frames/graphics of the animation. Let's say these are called p1, p2, ..., pn.  To animate, do:
plots:-display([p1, p2, ..., pn], insequence);

 

Solution...

Rouben Rostamian 7711
June 09 2018
2 0

I made some effort to understand the purpose of your calculations and untangle their logic. The result is what I have shown below.  Be sure to check every detail because it is possible that I have misinterpreted your worksheet or introduced typographical errors.

> 

restart;
> 

DEs := diff(x1(t), t) = x2(t),
       diff(x2(t), t) = l*u(x4(t))-k*sin(2*x1(t)),
       diff(x3(t), t) = 2*k*cos(2*x1(t))*x4(t),
       diff(x4(t), t) = -x3(t);

diff(x1(t), t) = x2(t), diff(x2(t), t) = l*u(x4(t))-k*sin(2*x1(t)), diff(x3(t), t) = 2*k*cos(2*x1(t))*x4(t), diff(x4(t), t) = -x3(t)

> 

u := x -> piecewise(x <= -1, -U_max, 1 < x, U_max, 0);

u := proc (x) options operator, arrow; piecewise(x <= -1, -U_max, 1 < x, U_max, 0) end proc

> 

bc := x1(0) = th0,
      x2(0) = q0,
      x3(t1) + 2*G1*x1(t1) = 0,
      x4(t1) + 2*G2*x2(t1) = 0;

x1(0) = th0, x2(0) = q0, x3(t1)+2*G1*x1(t1) = 0, x4(t1)+2*G2*x2(t1) = 0

Note that the parameters are specified as equations, not as assignments!

> 

params :=
  T0 = 300,
  G1 = 1.2e6,
  G2 = 1.2e6,
  k = 1.79333595e-6,
  l = 1/12e4,
  U_max = 5.0,
  th0 = .5,
  q0 = 0.1e-2;

T0 = 300, G1 = 0.12e7, G2 = 0.12e7, k = 0.179333595e-5, l = 0.8333333333e-5, U_max = 5.0, th0 = .5, q0 = 0.1e-2

> 

sys := eval({DEs,bc}, {params, t1=10.0}):    # note the t1=10.0
> 

dsol := dsolve(sys, numeric);

proc (x_bvp) local res, data, solnproc, _ndsol, outpoint, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then outpoint := evalf[_EnvDSNumericSaveDigits](x_bvp) else outpoint := evalf(x_bvp) end if; data := Array(1..4, {(1) = proc (outpoint) local X, Y, YP, yout, errproc, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; X := Vector(8, {(1) = .0, (2) = 1.428582867243206, (3) = 2.8572001869402115, (4) = 4.285793276873327, (5) = 5.714364424810432, (6) = 7.142918882256333, (7) = 8.571462217201669, (8) = 10.0}, datatype = float[8], order = C_order); Y := Matrix(8, 4, {(1, 1) = .5, (1, 2) = 0.10e-2, (1, 3) = -1218701.2042618706, (1, 4) = -12189157.44105805, (2, 1) = .5013845243918722, (2, 2) = 0.9383179909276191e-3, (2, 3) = -1218732.4725656176, (2, 4) = -10448118.857015757, (3, 1) = .5026809584542433, (3, 2) = 0.876630798947629e-3, (3, 3) = -1218758.8209411495, (3, 4) = -8706997.134761466, (4, 1) = .5038892428934152, (4, 2) = 0.8149412148685686e-3, (4, 3) = -1218780.2942966516, (4, 4) = -6965870.788812328, (5, 1) = .5050093791867353, (5, 2) = 0.753249394199044e-3, (5, 3) = -1218796.93510486, (5, 4) = -5224743.965523619, (6, 1) = .5060413706581044, (6, 2) = 0.6915553617706725e-3, (6, 3) = -1218808.7819558734, (6, 4) = -3483617.140681507, (7, 1) = .5069852198925933, (7, 2) = 0.6298591262419049e-3, (7, 3) = -1218815.8694750136, (7, 4) = -1742490.3515712114, (8, 1) = .5078409284568352, (8, 2) = 0.5681606938670946e-3, (8, 3) = -1218818.2282964045, (8, 4) = -1363.585665281027}, datatype = float[8], order = C_order); YP := Matrix(8, 4, {(1, 1) = 0.10e-2, (1, 2) = -0.431757068329379e-4, (1, 3) = -23.621210940059537, (1, 4) = 1218701.2042618706, (2, 1) = 0.9383179909276191e-3, (2, 2) = -0.4317838409608693e-4, (2, 3) = -20.15988007732689, (2, 4) = 1218732.4725656176, (3, 1) = 0.876630798947629e-3, (3, 2) = -0.43180880510326025e-4, (3, 3) = -16.732032844411012, (3, 4) = 1218758.8209411495, (4, 1) = 0.8149412148685686e-3, (4, 2) = -0.4318319801867119e-4, (4, 3) = -13.335132954895201, (4, 4) = 1218780.2942966516, (5, 1) = 0.753249394199044e-3, (5, 2) = -0.43185338547120616e-4, (5, 3) = -9.96647565352916, (5, 4) = 1218796.93510486, (6, 1) = 0.6915553617706725e-3, (6, 2) = -0.4318730388928721e-4, (6, 3) = -6.623331055162384, (6, 4) = 1218808.7819558734, (7, 1) = 0.6298591262419049e-3, (7, 2) = -0.4318909570042439e-4, (7, 3) = -3.302952234722308, (7, 4) = 1218815.8694750136, (8, 1) = 0.5681606938670946e-3, (8, 2) = -0.4319071549599433e-4, (8, 3) = -0.2577615512201638e-2, (8, 4) = 1218818.2282964045}, datatype = float[8], order = C_order); errproc := proc (x_bvp) local outpoint, X, Y, yout, L, V, i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; Digits := 15; outpoint := evalf(x_bvp); X := Vector(8, {(1) = .0, (2) = 1.428582867243206, (3) = 2.8572001869402115, (4) = 4.285793276873327, (5) = 5.714364424810432, (6) = 7.142918882256333, (7) = 8.571462217201669, (8) = 10.0}, datatype = float[8], order = C_order); Y := Matrix(8, 4, {(1, 1) = .0, (1, 2) = .0, (1, 3) = 0.39283088202469616e-9, (1, 4) = 0.1583150741983841e-8, (2, 1) = -0.7580062395781635e-16, (2, 2) = -0.3025210845716929e-18, (2, 3) = 0.37361311464173746e-9, (2, 4) = -0.17786774287792475e-8, (3, 1) = 0.2598274327370506e-15, (3, 2) = 0.33166261616014577e-18, (3, 3) = 0.474826688104923e-9, (3, 4) = 0.2021512657291743e-8, (4, 1) = -0.4601610101061015e-15, (4, 2) = -0.23276499915587446e-18, (4, 3) = 0.45038562915775664e-9, (4, 4) = 0.15693529661589289e-8, (5, 1) = -0.2840761139229594e-15, (5, 2) = 0.17448199333926153e-18, (5, 3) = -0.5548514105926573e-10, (5, 4) = 0.21764372136549502e-8, (6, 1) = -0.11898771241848955e-15, (6, 2) = -0.19825386038967814e-19, (6, 3) = 0.4353070721923952e-9, (6, 4) = -0.690459617210538e-9, (7, 1) = -0.23398097137467846e-15, (7, 2) = 0.2295525899574634e-18, (7, 3) = 0.4071210128732176e-9, (7, 4) = 0.5832347248695947e-9, (8, 1) = -0.12490890196393565e-15, (8, 2) = 0.18274259637436185e-18, (8, 3) = 0.6741001502551099e-9, (8, 4) = -0.3852601004217984e-12}, datatype = float[8], order = C_order); if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "right" then return X[8] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.1764372136549502e-9) elif outpoint = "errorproc" then error "this is already the error procedure" elif outpoint = "rawdata" then return [4, 8, [x1(t), x2(t), x3(t), x4(t)], X, Y] else return ('procname')(x_bvp) end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; V := array([1 = 4, 2 = 0]); if Digits <= trunc(evalhf(Digits)) then L := Vector(4, 'datatype' = 'float'[8]); yout := Vector(4, 'datatype' = 'float'[8]); evalhf(`dsolve/numeric/lagrange`(8, 4, X, Y, outpoint, var(yout), var(L), var(V))) else L := Vector(4, 'datatype' = 'sfloat'); yout := Vector(4, 'datatype' = 'sfloat'); `dsolve/numeric/lagrange`(8, 4, X, Y, outpoint, yout, L, V) end if; [t = outpoint, seq('[x1(t), x2(t), x3(t), x4(t)]'[i] = yout[i], i = 1 .. 4)] end proc; if not type(outpoint, 'numeric') then if outpoint = "start" or outpoint = "left" then return X[1] elif outpoint = "method" then return "bvp" elif outpoint = "right" then return X[8] elif outpoint = "order" then return 6 elif outpoint = "error" then return HFloat(2.1764372136549502e-9) elif outpoint = "errorproc" then return eval(errproc) elif outpoint = "rawdata" then return [4, 8, "depnames", X, Y, YP] else error "non-numeric value" end if end if; if outpoint < X[1] or X[8] < outpoint then error "solution is only defined in the range %1..%2", X[1], X[8] end if; if Digits <= trunc(evalhf(Digits)) and (_EnvInFsolve <> true or _EnvDSNumericSaveDigits <= trunc(evalhf(Digits))) then V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0, (3, 1) = .0, (3, 2) = .0, (4, 1) = .0, (4, 2) = .0, (5, 1) = .0, (5, 2) = .0, (6, 1) = .0, (6, 2) = .0, (7, 1) = .0, (7, 2) = .0}, datatype = float[8], order = C_order); yout := Vector(4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8]); evalhf(`dsolve/numeric/hermite`(8, 4, X, Y, YP, outpoint, var(yout), var(L), var(V))) else if _EnvInFsolve = true then Digits := _EnvDSNumericSaveDigits end if; V := array( 1 .. 6, [( 1 ) = (7), ( 2 ) = (0), ( 3 ) = (false), ( 4 ) = (false), ( 5 ) = (false), ( 6 ) = (false) ] ); L := Matrix(7, 2, {(1, 1) = 0., (1, 2) = 0., (2, 1) = 0., (2, 2) = 0., (3, 1) = 0., (3, 2) = 0., (4, 1) = 0., (4, 2) = 0., (5, 1) = 0., (5, 2) = 0., (6, 1) = 0., (6, 2) = 0., (7, 1) = 0., (7, 2) = 0.}, order = C_order); yout := Vector(4, {(1) = 0., (2) = 0., (3) = 0., (4) = 0.}); `dsolve/numeric/hermite`(8, 4, X, Y, YP, outpoint, yout, L, V) end if; [outpoint, seq(yout[i], i = 1 .. 4)] end proc, (2) = Array(0..0, {}), (3) = [t, x1(t), x2(t), x3(t), x4(t)], (4) = 0}); solnproc := data[1]; if not type(outpoint, 'numeric') then if outpoint = "solnprocedure" then return eval(solnproc) elif member(outpoint, ["start", "left", "right", "errorproc", "rawdata", "order", "error"]) then return solnproc(x_bvp) elif outpoint = "sysvars" then return data[3] elif procname <> unknown then return ('procname')(x_bvp) else _ndsol := pointto(data[2][0]); return ('_ndsol')(x_bvp) end if end if; try res := solnproc(outpoint); [t = res[1], seq('[x1(t), x2(t), x3(t), x4(t)]'[i] = res[i+1], i = 1 .. 4)] catch: error end try end proc

> 

dsol(0);

[t = 0., x1(t) = HFloat(0.4999999999999999), x2(t) = HFloat(9.999999999999998e-4), x3(t) = HFloat(-1218701.2042618704), x4(t) = HFloat(-1.2189157441058049e7)]

Conclusion:  At t=0 we have x3 = -1.21870120426187*10^6, x4 = -1.21891574410580*10^7

 

 

Download mw.mw

 

Newton's interpolation...

Rouben Rostamian 7711
June 05 2018
1 2

Polynomial interpolation through

Newton's forward difference scheme

> 

restart;

Procedure receives lists Xand Y of coordinates x__0, y__0, () .. (), x__k, y__k of points and
produces a polynomial of kth degree in the variable x (user-specified as the proc's
third argument) whose graph passes through those points.

The coordinates x__j must be distinct, but need not be equally spaced.

> 

N_interp := proc(X::list, Y::list, x)
  local i, j, k, N;
  k := nops(X);
  if nops(Y) <> k then
    error "The x and y arrays should be of equals lengths";
  end if;
  N := Matrix(k);
  for j from 1 to k do
    N[1,j] := Y[j];
  end do;
  for i from 2 to k do
    for j from 1 to k + 1 - i do
      N[i,j] := (N[i-1,j+1] - N[i-1,j])/(X[j+i-1]-X[j]);
    end do;
  end do;
  unassign('i', 'j');
  N[1,1] + add(N[i+1,1]*product(x-X[j], j=1..i), i=1..k-1);
end proc:

 

Example:

> 

X := [seq(i^2, i=0..3)];
Y := [seq(1/(i+1), i=0..3)];

[0, 1, 4, 9]

[1, 1/2, 1/3, 1/4]

> 

ans := N_interp(X,Y,lambda);   # asking for a polynomial in lambda

1-(1/2)*lambda+(1/9)*lambda*(lambda-1)-(17/1440)*lambda*(lambda-1)*(lambda-4)

> 

# verify the answer - should produce sequence of zeros if correct:
> 

seq(eval(ans, lambda=X[i]) - Y[i], i=1..nops(X));

0, 0, 0, 0

 

 

Download newton-interpolation.mw

Another way...

Rouben Rostamian 7711
June 04 2018
1 1

Tomleslie said there are other ways.  Here is one to produce the same object:

restart;
with(plots): with(plottools):
display(cuboid([0,0,0], [1,2,5]), scaling=constrained);

 

Obvsious solution...

Rouben Rostamian 7711
June 03 2018
0 4

You don't need Maple for that—just looking at the equations it is evident theta(r)≡0 and sigma(r)≡0 is a solution.

Solution...

Rouben Rostamian 7711
May 30 2018
1 1

I have fixed several issues with you worksheet as noted.  See if this does what you want.

> 

restart;
> 

C_t := -A(t)*omega[0]*(sigma-(diff(C(t), t)))-(1/4)*A(t)^5*beta[2]*omega[0]^2-(1/4)*A(t)^3*beta[1]*omega[0]^2-(1/2)*F[0]*cos(C(t))+5*alpha[2]*A(t)^5*(1/16)+3*alpha[1]*A(t)^3*(1/8) = 0;

-A(t)*omega[0]*(sigma-(diff(C(t), t)))-(1/4)*A(t)^5*beta[2]*omega[0]^2-(1/4)*A(t)^3*beta[1]*omega[0]^2-(1/2)*F[0]*cos(C(t))+(5/16)*alpha[2]*A(t)^5+(3/8)*alpha[1]*A(t)^3 = 0

 (1)
> 

A_t := (1/2)*mu*A(t)*omega[0]+(diff(A(t), t))*omega[0]-(1/2)*F[0]*sin(C(t)) = 0;

(1/2)*mu*A(t)*omega[0]+(diff(A(t), t))*omega[0]-(1/2)*F[0]*sin(C(t)) = 0

 (2)

We have a system of two first order differential equations, therefore we may
specify initial values, but not initial derivatives!

> 

ic := C(0)=0, A(0)=1;

C(0) = 0, A(0) = 1

 (3)

Note: params consists of equations, not assignments!

> 

params := F[0] =2,
          omega[0] = 1,            # note omega, not w !
          alpha[1] = 0.3331,
          alpha[2] = 0.1299,
          beta[1] = 0.3338,
          beta[2] = 0.1319,
          mu = 0.01,
          sigma = 0;

F[0] = 2, omega[0] = 1, alpha[1] = .3331, alpha[2] = .1299, beta[1] = .3338, beta[2] = .1319, mu = 0.1e-1, sigma = 0

 (4)
> 

sys := eval({C_t, A_t, ic}, {params});

{0.5000000000e-2*A(t)+diff(A(t), t)-sin(C(t)) = 0, A(t)*(diff(C(t), t))+0.761875000e-2*A(t)^5+0.4146250000e-1*A(t)^3-cos(C(t)) = 0, A(0) = 1, C(0) = 0}

 (5)
> 

dsol := dsolve(sys, numeric);

proc (x_rkf45) local _res, _dat, _vars, _solnproc, _xout, _ndsol, _pars, _n, _i; option `Copyright (c) 2000 by Waterloo Maple Inc. All rights reserved.`; if 1 < nargs then error "invalid input: too many arguments" end if; _EnvDSNumericSaveDigits := Digits; Digits := 15; if _EnvInFsolve = true then _xout := evalf[_EnvDSNumericSaveDigits](x_rkf45) else _xout := evalf(x_rkf45) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _dtbl, _dat, _vmap, _x0, _y0, _val, _dig, _n, _ne, _nd, _nv, _pars, _ini, _par, _i, _j, _k, _src; option `Copyright (c) 2002 by Waterloo Maple Inc. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _pars := []; _dtbl := array( 1 .. 4, [( 1 ) = (array( 1 .. 26, [( 1 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 2 ) = (datatype = float[8], order = C_order, storage = rectangular), ( 3 ) = ([0, 0, 0, Array(1..0, 1..2, {}, datatype = float[8], order = C_order)]), ( 4 ) = (Array(1..63, {(1) = 2, (2) = 2, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 1, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 1, (19) = 30000, (20) = 0, (21) = 0, (22) = 1, (23) = 4, (24) = 0, (25) = 1, (26) = 15, (27) = 1, (28) = 0, (29) = 1, (30) = 3, (31) = 3, (32) = 0, (33) = 1, (34) = 0, (35) = 0, (36) = 0, (37) = 0, (38) = 0, (39) = 0, (40) = 0, (41) = 0, (42) = 0, (43) = 1, (44) = 0, (45) = 0, (46) = 0, (47) = 0, (48) = 0, (49) = 0, (50) = 50, (51) = 1, (52) = 0, (53) = 0, (54) = 0, (55) = 0, (56) = 0, (57) = 0, (58) = 0, (59) = 10000, (60) = 0, (61) = 1000, (62) = 0, (63) = 0}, datatype = integer[8])), ( 5 ) = (Array(1..28, {(1) = .0, (2) = 0.10e-5, (3) = .0, (4) = 0.500001e-14, (5) = .0, (6) = 0.5308191426282788e-2, (7) = .0, (8) = 0.10e-5, (9) = .0, (10) = .0, (11) = .0, (12) = .0, (13) = 1.0, (14) = .0, (15) = .49999999999999, (16) = .0, (17) = 1.0, (18) = 1.0, (19) = .0, (20) = .0, (21) = 1.0, (22) = 1.0, (23) = .0, (24) = .0, (25) = 0.10e-14, (26) = .0, (27) = .0, (28) = .0}, datatype = float[8], order = C_order)), ( 6 ) = (Array(1..2, {(1) = 1.0, (2) = .0}, datatype = float[8], order = C_order)), ( 7 ) = ([Array(1..4, 1..7, {(1, 1) = .0, (1, 2) = .203125, (1, 3) = .3046875, (1, 4) = .75, (1, 5) = .8125, (1, 6) = .40625, (1, 7) = .8125, (2, 1) = 0.6378173828125e-1, (2, 2) = .0, (2, 3) = .279296875, (2, 4) = .27237892150878906, (2, 5) = -0.9686851501464844e-1, (2, 6) = 0.1956939697265625e-1, (2, 7) = .5381584167480469, (3, 1) = 0.31890869140625e-1, (3, 2) = .0, (3, 3) = -.34375, (3, 4) = -.335235595703125, (3, 5) = .2296142578125, (3, 6) = .41748046875, (3, 7) = 11.480712890625, (4, 1) = 0.9710520505905151e-1, (4, 2) = .0, (4, 3) = .40350341796875, (4, 4) = 0.20297467708587646e-1, (4, 5) = -0.6054282188415527e-2, (4, 6) = -0.4770040512084961e-1, (4, 7) = .77858567237854}, datatype = float[8], order = C_order), Array(1..6, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = 1.0, (2, 1) = .25, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = 1.0, (3, 1) = .1875, (3, 2) = .5625, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (3, 6) = 2.0, (4, 1) = .23583984375, (4, 2) = -.87890625, (4, 3) = .890625, (4, 4) = .0, (4, 5) = .0, (4, 6) = .2681884765625, (5, 1) = .1272735595703125, (5, 2) = -.5009765625, (5, 3) = .44921875, (5, 4) = -0.128936767578125e-1, (5, 5) = .0, (5, 6) = 0.626220703125e-1, (6, 1) = -0.927734375e-1, (6, 2) = .626220703125, (6, 3) = -.4326171875, (6, 4) = .1418304443359375, (6, 5) = -0.861053466796875e-1, (6, 6) = .3131103515625}, datatype = float[8], order = C_order), Array(1..6, {(1) = .0, (2) = .386, (3) = .21, (4) = .63, (5) = 1.0, (6) = 1.0}, datatype = float[8], order = C_order), Array(1..6, {(1) = .25, (2) = -.1043, (3) = .1035, (4) = -0.362e-1, (5) = .0, (6) = .0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 1.544, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = .9466785280815533, (3, 2) = .25570116989825814, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = 3.3148251870684886, (4, 2) = 2.896124015972123, (4, 3) = .9986419139977808, (4, 4) = .0, (4, 5) = .0, (5, 1) = 1.2212245092262748, (5, 2) = 6.019134481287752, (5, 3) = 12.537083329320874, (5, 4) = -.687886036105895, (5, 5) = .0, (6, 1) = 1.2212245092262748, (6, 2) = 6.019134481287752, (6, 3) = 12.537083329320874, (6, 4) = -.687886036105895, (6, 5) = 1.0}, datatype = float[8], order = C_order), Array(1..6, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = -5.6688, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (3, 1) = -2.4300933568337584, (3, 2) = -.20635991570891224, (3, 3) = .0, (3, 4) = .0, (3, 5) = .0, (4, 1) = -.10735290581452621, (4, 2) = -9.594562251021896, (4, 3) = -20.470286148096154, (4, 4) = .0, (4, 5) = .0, (5, 1) = 7.496443313968615, (5, 2) = -10.246804314641219, (5, 3) = -33.99990352819906, (5, 4) = 11.708908932061595, (5, 5) = .0, (6, 1) = 8.083246795922411, (6, 2) = -7.981132988062785, (6, 3) = -31.52159432874373, (6, 4) = 16.319305431231363, (6, 5) = -6.0588182388340535}, datatype = float[8], order = C_order), Array(1..3, 1..5, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (2, 1) = 10.126235083446911, (2, 2) = -7.487995877607633, (2, 3) = -34.800918615557414, (2, 4) = -7.9927717075687275, (2, 5) = 1.0251377232956207, (3, 1) = -.6762803392806898, (3, 2) = 6.087714651678606, (3, 3) = 16.43084320892463, (3, 4) = 24.767225114183653, (3, 5) = -6.5943891257167815}, datatype = float[8], order = C_order)]), ( 9 ) = ([Array(1..2, {(1) = .1, (2) = .1}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..2, {(1, 1) = .0, (1, 2) = .0, (2, 1) = .0, (2, 2) = .0}, datatype = float[8], order = C_order), Array(1..2, 1..6, {(1, 1) = .0, (1, 2) = .0, (1, 3) = .0, (1, 4) = .0, (1, 5) = .0, (1, 6) = .0, (2, 1) = .0, (2, 2) = .0, (2, 3) = .0, (2, 4) = .0, (2, 5) = .0, (2, 6) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = 0, (2) = 0}, datatype = integer[8]), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..4, {(1) = .0, (2) = .0, (3) = .0, (4) = .0}, datatype = float[8], order = C_order)]), ( 8 ) = ([Array(1..2, {(1) = 1.0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = .0, (2) = .0}, datatype = float[8], order = C_order), Array(1..2, {(1) = -0.50e-2, (2) = .95091875}, datatype = float[8], order = C_order), 0, 0]), ( 11 ) = (Array(1..6, 0..2, {(1, 1) = .0, (1, 2) = .0, (2, 0) = .0, (2, 1) = .0, (2, 2) = .0, (3, 0) = .0, (3, 1) = .0, (3, 2) = .0, (4, 0) = .0, (4, 1) = .0, (4, 2) = .0, (5, 0) = .0, (5, 1) = .0, (5, 2) = .0, (6, 0) = .0, (6, 1) = .0, (6, 2) = .0}, datatype = float[8], order = C_order)), ( 10 ) = ([proc (N, X, Y, YP) option `[Y[1] = A(t), Y[2] = C(t)]`; YP[1] := -0.5000000000e-2*Y[1]+sin(Y[2]); YP[2] := (-0.761875000e-2*Y[1]^5-0.4146250000e-1*Y[1]^3+cos(Y[2]))/Y[1]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 13 ) = (), ( 12 ) = (), ( 15 ) = ("rkf45"), ( 14 ) = ([0, 0]), ( 18 ) = ([]), ( 19 ) = (0), ( 16 ) = ([0, 0, 0, []]), ( 17 ) = ([proc (N, X, Y, YP) option `[Y[1] = A(t), Y[2] = C(t)]`; YP[1] := -0.5000000000e-2*Y[1]+sin(Y[2]); YP[2] := (-0.761875000e-2*Y[1]^5-0.4146250000e-1*Y[1]^3+cos(Y[2]))/Y[1]; 0 end proc, -1, 0, 0, 0, 0, 0, 0]), ( 22 ) = (0), ( 23 ) = (0), ( 20 ) = ([]), ( 21 ) = (0), ( 26 ) = (Array(1..0, {})), ( 25 ) = (Array(1..0, {})), ( 24 ) = (0) ] )) ] ); _y0 := Array(0..2, {(1) = 0., (2) = 1.}); _vmap := array( 1 .. 2, [( 1 ) = (1), ( 2 ) = (2) ] ); _x0 := _dtbl[1][5][5]; _n := _dtbl[1][4][1]; _ne := _dtbl[1][4][3]; _nd := _dtbl[1][4][4]; _nv := _dtbl[1][4][16]; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then if _Env_smart_dsolve_numeric = true or _dtbl[1][4][10] = 1 then if _xout = "left" then if type(_dtbl[2], 'table') then return _dtbl[2][5][1] end if elif _xout = "right" then if type(_dtbl[3], 'table') then return _dtbl[3][5][1] end if end if end if; return _dtbl[1][5][5] elif _xout = "method" then return _dtbl[1][15] elif _xout = "storage" then return evalb(_dtbl[1][4][10] = 1) elif _xout = "leftdata" then if not type(_dtbl[2], 'array') then return NULL else return eval(_dtbl[2]) end if elif _xout = "rightdata" then if not type(_dtbl[3], 'array') then return NULL else return eval(_dtbl[3]) end if elif _xout = "enginedata" then return eval(_dtbl[1]) elif _xout = "enginereset" then _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); return NULL elif _xout = "initial" then return procname(_y0[0]) elif _xout = "laxtol" then return _dtbl[`if`(member(_dtbl[4], {2, 3}), _dtbl[4], 1)][5][18] elif _xout = "numfun" then return `if`(member(_dtbl[4], {2, 3}), _dtbl[_dtbl[4]][4][18], 0) elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return procname(_y0[0]), [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _dtbl[4] <> 2 and _dtbl[4] <> 3 or _x0-_dtbl[_dtbl[4]][5][1] = 0. then error "no information is available on last computed point" else _xout := _dtbl[_dtbl[4]][5][1] end if elif _xout = "function" then if _dtbl[1][4][33]-2. = 0 then return eval(_dtbl[1][10], 1) else return eval(_dtbl[1][10][1], 1) end if elif _xout = "map" then return copy(_vmap) elif type(_xin, `=`) and type(rhs(_xin), 'list') and member(lhs(_xin), {"initial", "parameters", "initial_and_parameters"}) then _ini, _par := [], []; if lhs(_xin) = "initial" then _ini := rhs(_xin) elif lhs(_xin) = "parameters" then _par := rhs(_xin) elif select(type, rhs(_xin), `=`) <> [] then _par, _ini := selectremove(type, rhs(_xin), `=`) elif nops(rhs(_xin)) < nops(_pars)+1 then error "insufficient data for specification of initial and parameters" else _par := rhs(_xin)[-nops(_pars) .. -1]; _ini := rhs(_xin)[1 .. -nops(_pars)-1] end if; _xout := lhs(_xout); if _par <> [] then `dsolve/numeric/process_parameters`(_n, _pars, _par, _y0) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n-_ne, _ini, _y0, _pars, _vmap) end if; `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars); if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') and _dtbl[1][4][10] <> 1 then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)] elif _xout = "parameters" then return [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] else return [_y0[0], seq(_y0[_vmap[_i]], _i = 1 .. _n-_ne)], [seq(_y0[_n+_i], _i = 1 .. nops(_pars))] end if elif _xin = "eventstop" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then return 0 end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 <= _dtbl[5-_i][4][9] then _i := 5-_i; _dtbl[4] := _i; _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) elif 100 <= _dtbl[_i][4][9] then _j := round(_dtbl[_i][4][17]); return round(_dtbl[_i][3][1][_j, 1]) else return 0 end if elif _xin = "eventstatus" then if _nv = 0 then error "this solution has no events" end if; _i := [selectremove(proc (a) options operator, arrow; _dtbl[1][3][1][a, 7] = 1 end proc, {seq(_j, _j = 1 .. round(_dtbl[1][3][1][_nv+1, 1]))})]; return ':-enabled' = _i[1], ':-disabled' = _i[2] elif _xin = "eventclear" then if _nv = 0 then error "this solution has no events" end if; _i := _dtbl[4]; if _i <> 2 and _i <> 3 then error "no events to clear" end if; if _dtbl[_i][4][10] = 1 and assigned(_dtbl[5-_i]) and _dtbl[_i][4][9] < 100 and 100 < _dtbl[5-_i][4][9] then _dtbl[4] := 5-_i; _i := 5-_i end if; if _dtbl[_i][4][9] < 100 then error "no events to clear" elif _nv < _dtbl[_i][4][9]-100 then error "event error condition cannot be cleared" else _j := _dtbl[_i][4][9]-100; if irem(round(_dtbl[_i][3][1][_j, 4]), 2) = 1 then error "retriggerable events cannot be cleared" end if; _j := round(_dtbl[_i][3][1][_j, 1]); for _k to _nv do if _dtbl[_i][3][1][_k, 1] = _j then if _dtbl[_i][3][1][_k, 2] = 3 then error "range events cannot be cleared" end if; _dtbl[_i][3][1][_k, 8] := _dtbl[_i][3][1][_nv+1, 8] end if end do; _dtbl[_i][4][17] := 0; _dtbl[_i][4][9] := 0; if _dtbl[1][4][10] = 1 then if _i = 2 then try procname(procname("left")) catch: end try else try procname(procname("right")) catch: end try end if end if end if; return elif type(_xin, `=`) and member(lhs(_xin), {"eventdisable", "eventenable"}) then if _nv = 0 then error "this solution has no events" end if; if type(rhs(_xin), {('list')('posint'), ('set')('posint')}) then _i := {op(rhs(_xin))} elif type(rhs(_xin), 'posint') then _i := {rhs(_xin)} else error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; if select(proc (a) options operator, arrow; _nv < a end proc, _i) <> {} then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _k := {}; for _j to _nv do if member(round(_dtbl[1][3][1][_j, 1]), _i) then _k := `union`(_k, {_j}) end if end do; _i := _k; if lhs(_xin) = "eventdisable" then _dtbl[4] := 0; _j := [evalb(assigned(_dtbl[2]) and member(_dtbl[2][4][17], _i)), evalb(assigned(_dtbl[3]) and member(_dtbl[3][4][17], _i))]; for _k in _i do _dtbl[1][3][1][_k, 7] := 0; if assigned(_dtbl[2]) then _dtbl[2][3][1][_k, 7] := 0 end if; if assigned(_dtbl[3]) then _dtbl[3][3][1][_k, 7] := 0 end if end do; if _j[1] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[2][3][4][_k, 1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to defined init `, _dtbl[2][3][4][_k, 1]); _dtbl[2][3][1][_k, 8] := _dtbl[2][3][4][_k, 1] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to rate hysteresis init `, _dtbl[2][5][24]); _dtbl[2][3][1][_k, 8] := _dtbl[2][5][24] elif _dtbl[2][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[2][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to initial init `, _x0); _dtbl[2][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #2, event code `, _k, ` to fireinitial init `, _x0-1); _dtbl[2][3][1][_k, 8] := _x0-1 end if end do; _dtbl[2][4][17] := 0; _dtbl[2][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("left")) end if end if; if _j[2] then for _k to _nv+1 do if _k <= _nv and not type(_dtbl[3][3][4][_k, 2], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to defined init `, _dtbl[3][3][4][_k, 2]); _dtbl[3][3][1][_k, 8] := _dtbl[3][3][4][_k, 2] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to rate hysteresis init `, _dtbl[3][5][24]); _dtbl[3][3][1][_k, 8] := _dtbl[3][5][24] elif _dtbl[3][3][1][_k, 2] = 0 and irem(iquo(round(_dtbl[3][3][1][_k, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to initial init `, _x0); _dtbl[3][3][1][_k, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #3, event code `, _k, ` to fireinitial init `, _x0+1); _dtbl[3][3][1][_k, 8] := _x0+1 end if end do; _dtbl[3][4][17] := 0; _dtbl[3][4][9] := 0; if _dtbl[1][4][10] = 1 then procname(procname("right")) end if end if else for _k in _i do _dtbl[1][3][1][_k, 7] := 1 end do; _dtbl[2] := evaln(_dtbl[2]); _dtbl[3] := evaln(_dtbl[3]); _dtbl[4] := 0; if _dtbl[1][4][10] = 1 then if _x0 <= procname("right") then try procname(procname("right")) catch: end try end if; if procname("left") <= _x0 then try procname(procname("left")) catch: end try end if end if end if; return elif type(_xin, `=`) and lhs(_xin) = "eventfired" then if not type(rhs(_xin), 'list') then error "'eventfired' must be specified as a list" end if; if _nv = 0 then error "this solution has no events" end if; if _dtbl[4] <> 2 and _dtbl[4] <> 3 then error "'direction' must be set prior to calling/setting 'eventfired'" end if; _i := _dtbl[4]; _val := NULL; if not assigned(_EnvEventRetriggerWarned) then _EnvEventRetriggerWarned := false end if; for _k in rhs(_xin) do if type(_k, 'integer') then _src := _k elif type(_k, 'integer' = 'anything') and type(evalf(rhs(_k)), 'numeric') then _k := lhs(_k) = evalf[max(Digits, 18)](rhs(_k)); _src := lhs(_k) else error "'eventfired' entry is not valid: %1", _k end if; if _src < 1 or round(_dtbl[1][3][1][_nv+1, 1]) < _src then error "event identifiers must be integers in the range 1..%1", round(_dtbl[1][3][1][_nv+1, 1]) end if; _src := {seq(`if`(_dtbl[1][3][1][_j, 1]-_src = 0., _j, NULL), _j = 1 .. _nv)}; if nops(_src) <> 1 then error "'eventfired' can only be set/queried for root-finding events and time/interval events" end if; _src := _src[1]; if _dtbl[1][3][1][_src, 2] <> 0. and _dtbl[1][3][1][_src, 2]-2. <> 0. then error "'eventfired' can only be set/queried for root-finding events and time/interval events" elif irem(round(_dtbl[1][3][1][_src, 4]), 2) = 1 then if _EnvEventRetriggerWarned = false then WARNING(`'eventfired' has no effect on events that retrigger`) end if; _EnvEventRetriggerWarned := true end if; if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then _val := _val, undefined elif type(_dtbl[_i][3][4][_src, _i-1], 'undefined') or _i = 2 and _dtbl[2][3][1][_src, 8] < _dtbl[2][3][4][_src, 1] or _i = 3 and _dtbl[3][3][4][_src, 2] < _dtbl[3][3][1][_src, 8] then _val := _val, _dtbl[_i][3][1][_src, 8] else _val := _val, _dtbl[_i][3][4][_src, _i-1] end if; if type(_k, `=`) then if _dtbl[_i][3][1][_src, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_src, 4]), 32), 2) = 1 then error "cannot set event code for a rate hysteresis event" end if; userinfo(3, {'events', 'eventreset'}, `manual set event code `, _src, ` to value `, rhs(_k)); _dtbl[_i][3][1][_src, 8] := rhs(_k); _dtbl[_i][3][4][_src, _i-1] := rhs(_k) end if end do; return [_val] elif type(_xin, `=`) and lhs(_xin) = "direction" then if not member(rhs(_xin), {-1, 1, ':-left', ':-right'}) then error "'direction' must be specified as either '1' or 'right' (positive) or '-1' or 'left' (negative)" end if; _src := `if`(_dtbl[4] = 2, -1, `if`(_dtbl[4] = 3, 1, undefined)); _i := `if`(member(rhs(_xin), {1, ':-right'}), 3, 2); _dtbl[4] := _i; _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #4, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if; return _src elif _xin = "eventcount" then if _dtbl[1][3][1] = 0 or _dtbl[4] <> 2 and _dtbl[4] <> 3 then return 0 else return round(_dtbl[_dtbl[4]][3][1][_nv+1, 12]) end if else return "procname" end if end if; if _xout = _x0 then return [_x0, seq(evalf(_dtbl[1][6][_vmap[_i]]), _i = 1 .. _n-_ne)] end if; _i := `if`(_x0 <= _xout, 3, 2); if _xin = "last" and 0 < _dtbl[_i][4][9] and _dtbl[_i][4][9] < 100 then _dat := eval(_dtbl[_i], 2); _j := _dat[4][20]; return [_dat[11][_j, 0], seq(_dat[11][_j, _vmap[_i]], _i = 1 .. _n-_ne-_nd), seq(_dat[8][1][_vmap[_i]], _i = _n-_ne-_nd+1 .. _n-_ne)] end if; if not type(_dtbl[_i], 'array') then _dtbl[_i] := `dsolve/numeric/SC/IVPdcopy`(_dtbl[1], `if`(assigned(_dtbl[_i]), _dtbl[_i], NULL)); if 0 < _nv then for _j to _nv+1 do if _j <= _nv and not type(_dtbl[_i][3][4][_j, _i-1], 'undefined') then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to defined init `, _dtbl[_i][3][4][_j, _i-1]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][3][4][_j, _i-1] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 32), 2) = 1 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to rate hysteresis init `, _dtbl[_i][5][24]); _dtbl[_i][3][1][_j, 8] := _dtbl[_i][5][24] elif _dtbl[_i][3][1][_j, 2] = 0 and irem(iquo(round(_dtbl[_i][3][1][_j, 4]), 2), 2) = 0 then userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to initial init `, _x0); _dtbl[_i][3][1][_j, 8] := _x0 else userinfo(3, {'events', 'eventreset'}, `reinit #5, event code `, _j, ` to fireinitial init `, _x0-2*_i+5.0); _dtbl[_i][3][1][_j, 8] := _x0-2*_i+5.0 end if end do end if end if; if _xin <> "last" then if 0 < 0 then if `dsolve/numeric/checkglobals`(op(_dtbl[1][14]), _pars, _n, _y0) then `dsolve/numeric/SC/reinitialize`(_dtbl, _y0, _n, procname, _pars, _i) end if end if; if _dtbl[1][4][7] = 0 then error "parameters must be initialized before solution can be computed" end if end if; _dat := eval(_dtbl[_i], 2); _dtbl[4] := _i; try _src := `dsolve/numeric/SC/IVPrun`(_dat, _xout) catch: userinfo(2, `dsolve/debug`, print(`Exception in solnproc:`, [lastexception][2 .. -1])); error end try; if _dat[17] <> _dtbl[1][17] then _dtbl[1][17] := _dat[17]; _dtbl[1][10] := _dat[10] end if; if _src = 0 and 100 < _dat[4][9] then _val := _dat[3][1][_nv+1, 8] else _val := _dat[11][_dat[4][20], 0] end if; if _src <> 0 or _dat[4][9] <= 0 then _dtbl[1][5][1] := _xout else _dtbl[1][5][1] := _val end if; if _i = 3 and _val < _xout then Rounding := -infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further right of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further right of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further right of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further right of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further right of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further right of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further right of %1", evalf[8](_val) end if elif _i = 2 and _xout < _val then Rounding := infinity; if _dat[4][9] = 1 then error "cannot evaluate the solution further left of %1, probably a singularity", evalf[8](_val) elif _dat[4][9] = 2 then error "cannot evaluate the solution further left of %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_val) elif _dat[4][9] = 3 then if _dat[4][25] = 3 then error "cannot evaluate the solution past the initial point, problem may be initially singular or improperly set up" else error "cannot evaluate the solution past the initial point, problem may be complex, initially singular or improperly set up" end if elif _dat[4][9] = 4 then error "cannot evaluate the solution further left of %1, accuracy goal cannot be achieved with specified 'minstep'", evalf[8](_val) elif _dat[4][9] = 5 then error "cannot evaluate the solution further left of %1, too many step failures, tolerances may be too loose for problem", evalf[8](_val) elif _dat[4][9] = 6 then error "cannot evaluate the solution further left of %1, cannot downgrade delay storage for problems with delay derivative order > 1, try increasing delaypts", evalf[8](_val) elif _dat[4][9] = 10 then error "cannot evaluate the solution further right of %1, interrupt requested", evalf[8](_val) elif 100 < _dat[4][9] then if _dat[4][9]-100 = _nv+1 then error "constraint projection failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+2 then error "index-1 and derivative evaluation failure on event at t=%1", evalf[8](_val) elif _dat[4][9]-100 = _nv+3 then error "maximum number of event iterations reached (%1) at t=%2", round(_dat[3][1][_nv+1, 3]), evalf[8](_val) else if _Env_dsolve_nowarnstop <> true then `dsolve/numeric/warning`(StringTools:-FormatMessage("cannot evaluate the solution further left of %1, event #%2 triggered a halt", evalf[8](_val), round(_dat[3][1][_dat[4][9]-100, 1]))) end if; Rounding := 'nearest'; _xout := _val end if else error "cannot evaluate the solution further left of %1", evalf[8](_val) end if end if; if _EnvInFsolve = true then _dig := _dat[4][26]; _dat[4][26] := _EnvDSNumericSaveDigits; _Env_dsolve_SC_native := true; if _dat[4][25] = 1 then _i := 1; _dat[4][25] := 2 else _i := _dat[4][25] end if; _val := `dsolve/numeric/SC/IVPval`(_dat, _xout, _src); _dat[4][25] := _i; _dat[4][26] := _dig; [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] else Digits := _dat[4][26]; _val := `dsolve/numeric/SC/IVPval`(eval(_dat, 2), _xout, _src); [_xout, seq(_val[_vmap[_i]], _i = 1 .. _n-_ne)] end if end proc, (2) = Array(0..0, {}), (3) = [t, A(t), C(t)], (4) = []}); _vars := _dat[3]; _pars := map(rhs, _dat[4]); _n := nops(_vars)-1; _solnproc := _dat[1]; if not type(_xout, 'numeric') then if member(x_rkf45, ["start", 'start', "method", 'method', "left", 'left', "right", 'right', "leftdata", "rightdata", "enginedata", "eventstop", 'eventstop', "eventclear", 'eventclear', "eventstatus", 'eventstatus', "eventcount", 'eventcount', "laxtol", 'laxtol', "numfun", 'numfun', NULL]) then _res := _solnproc(convert(x_rkf45, 'string')); if 1 < nops([_res]) then return _res elif type(_res, 'array') then return eval(_res, 1) elif _res <> "procname" then return _res end if elif member(x_rkf45, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_rkf45, 'string'); _res := _solnproc(_xout); if _xout = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] end if elif type(_xout, `=`) and member(lhs(_xout), ["initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(lhs(x_rkf45), 'string') = rhs(x_rkf45); if type(rhs(_xout), 'list') then _res := _solnproc(_xout) else error "initial and/or parameter values must be specified in a list" end if; if lhs(_xout) = "initial" then return [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] elif lhs(_xout) = "parameters" then return [seq(_pars[_i] = _res[_i], _i = 1 .. nops(_pars))] else return [seq(_vars[_i+1] = [_res][1][_i+1], _i = 0 .. _n), seq(_pars[_i] = [_res][2][_i], _i = 1 .. nops(_pars))] end if elif type(_xout, `=`) and member(lhs(_xout), ["eventdisable", 'eventdisable', "eventenable", 'eventenable', "eventfired", 'eventfired', "direction", 'direction', NULL]) then return _solnproc(convert(lhs(x_rkf45), 'string') = rhs(x_rkf45)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_rkf45) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_rkf45) end if end if; try _res := _solnproc(_xout); [seq(_vars[_i+1] = _res[_i+1], _i = 0 .. _n)] catch: error end try end proc

 (6)
> 

plots:-odeplot(dsol, [t, A(t)], t=0..20);

 
> 

plots:-odeplot(dsol, [t, C(t)], t=0..20);

 
> 

plots:-odeplot(dsol, [C(t), A(t)], t=0..20);

 
 

 

Download mw.mw

Use Explore...

Rouben Rostamian 7711
May 25 2018
0 3

Maple's Explore is the right tool for solving your problem.  Explore assigns graphical "sliders" to the problem's parameters.  You may change the parameter values through the sliders and observe the effect on the problem's solutions.

Your worksheet contains seven equations, while you said you are interested in solving three of them,  I have extracted those three equations into a new worksheet and have added the Explore command.  I have assumed that the parameters take values in the range -5 to +5.  You may change those ranges as needed.

Here is the worksheet:  explore.mw

 

Do it like this...

Rouben Rostamian 7711
May 24 2018
1 0
> 

restart;
> 

f := sin(x);

sin(x)

> 

solve(f,x);

0

> 

solve(f,x,allsolutions);

Pi*_Z1

> 

about(_Z1);

Originally _Z1, renamed _Z1~:
  is assumed to be: integer

Solution...

Rouben Rostamian 7711
May 23 2018
1 1

In your worksheet you are attempting to solve a system of six differential equations by splitting the process into two stages.  In the first stage you solve for the unknowns p, q, r, and then in the second stage you solve for theta, phi and psi which depend on p, q, r.

The preferred way is to solve the entire system of six equations in one pass.  I have modified your worksheet to do that.  Along the way I have made several other minor adjustments.  I believe these should be acceptable to your Maple 13, Let us know if it doesn't work.

Пример02-новый.mw

 

Command-line tools...

Rouben Rostamian 7711
May 22 2018
4 0

I make a lot of animations with varying degrees of complexity.  I use Maple to generate frames, and I use ImageMagick to assemble the frames into an animation with complete control on everything.

ImageMagick is an extensive suite of software tools for image manipulation.  It is open and free.  I use it on Linux, but I understand that it also works on other platforms.  In all cases it is accessed through the command-line.

Here I will describe the steps needed to accomplish what you have asked.

  1. Save your Maple animation to an (animated) GIF file as you have described.  Let's call that file anim-raw.gif.
  2. Apply ImageMagick's convert command to split that file into individual frames:
          convert   anim-raw.gif   frame-%03d.gif
    This will write the individual frames into files named frame-000.gif through frame-nnn.gif.  If you have more that a 1000 frames, change the %03d to %04d.  That will write the files frame-0000.gif, etc.
  3. Now put together the individual frames into a single animated GIF called anim.gif:
          convert   -delay 10   -loop 2   frame-*.gif   anim.gif
    The "-loop 2" option says that we want to play the animation twice and then stop.  Change the 2 to anything you wish, or say "-loop 0" if you want the animation to play forever.
    The "-delay 10" option specifies the delay between frames.  The delay is measured in units of 1/100 of a second, so "-delay 10" means 1/10 of a second between frames.

I think that addresses all the questions that you have asked, and as we see, it takes just two commands typed on the command-line.

There are many more options available.  Consult ImageMagick's documentation.

 

 

 

 

See if this is good...

Rouben Rostamian 7711
May 21 2018
1 1

It seems to me that your problem originates in TeXShop.  I don't have a Mac and I have no experience with TeXShop.  I do all my work on Linux, and with TeXLive which comes with it.

I downloaded your Konstruktion_Masterfile.mws, loaded it into Maple 2018, and saved the result as LaTeX through File->Export As...

Then compiled the resulting file through the command-line:

pdflatex Konstruktion_Masterfile.tex      # first pass
pdflatex Konstruktion_Masterfile.tex      # second pass

There were 6 error messages which I ingored.  They all were about missing { and } on line 4779.  I didn't try to fix those; you should be able to do it since you are familiar with the content. 

Here are Konstruktion_Masterfile.tex (which I have renamed Konstruktion_Masterfile.txt) and Konstruktion_Masterfile.pdf.  They look OK to me.

Konstruktion_Masterfile.txt
Konstruktion_Masterfile.pdf

 

Solution...

Rouben Rostamian 7711
May 10 2018
1 1
> 

restart;

I will solve the PDE in the domain 0 < x and x < 1, 0 < y and y < 1.  Solving in a general rectangular domain is no more difficult.

> 

pde := diff(u(x, y),x,x) + diff(u(x,y),y,y) = 2*x^3 + 6*x*y*(1-y);

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 2*x^3+6*x*y*(1-y)

> 

bc1 := u(0,y) = 0, u(1, y) = y*(1-y);

u(0, y) = 0, u(1, y) = y*(1-y)

> 

bc2 := u(x,0) = 0, u(x,1) = 0;

u(x, 0) = 0, u(x, 1) = 0

In view of bc2, the solution may be expressed in a Fourier series in the y variable, where the coefficients f__k(x) are to be determined:

> 

sol := u(x,y) = Sum(f[k](x)*sin(k*Pi*y), k=1..infinity);

u(x, y) = Sum(f[k](x)*sin(k*Pi*y), k = 1 .. infinity)

The boundary conditions bc2 are already satisfied.  This solution candidate should satisfy the pde and the boundary conditions bc1.

Step 1:  Let's look at bc1.  The first equation in bc1 says that:

> 

eval(rhs(sol), x=0) = rhs(bc1[1]);

Sum(f[k](0)*sin(k*Pi*y), k = 1 .. infinity) = 0

which then implies that

> 

ic1 := f[k](0) = 0;

f[k](0) = 0

Step 2:  The second equation in bc1 says that

> 

tmp := eval(rhs(sol), x=1) = rhs(bc1[2]);

Sum(f[k](1)*sin(k*Pi*y), k = 1 .. infinity) = y*(1-y)

To determine the f__k(1) from here, we multiply the equation by sin(n*Pi*y) and integrate:

> 

Int(lhs(tmp)*sin(n*Pi*y), y=0..1) = Int(rhs(tmp)*sin(n*Pi*y), y=0..1);

Int((Sum(f[k](1)*sin(k*Pi*y), k = 1 .. infinity))*sin(n*Pi*y), y = 0 .. 1) = Int(y*(1-y)*sin(n*Pi*y), y = 0 .. 1)

The sum of the left-hand side reduces to a single term because the integral of sin(k*Pi*y)*sin(n*Pi*y) is zero when k <> n.  The integral on the right-hand side is easy to calculate.  We thus arrive at

> 

int(f[k](1)*sin(k*Pi*y)^2, y=0..1) = int(y*(1-y)*sin(k*Pi*y), y=0..1) assuming k::integer:
ic2 := isolate(%, f[k](1));

f[k](1) = 4*(-(-1)^k+1)/(k^3*Pi^3)

Step 3:  Plug the solution candidate u(x, y) into he pde:

> 

tmp := eval(pde, sol);

Sum((diff(diff(f[k](x), x), x))*sin(k*Pi*y), k = 1 .. infinity)+Sum(-f[k](x)*k^2*Pi^2*sin(k*Pi*y), k = 1 .. infinity) = 2*x^3+6*x*y*(1-y)

Multiply this by sin(n*Pi*y) and integrate.

> 

Int(lhs(tmp)*sin(n*Pi*y), y=0..1) = Int(rhs(tmp)*sin(n*Pi*y), y=0..1);

Int((Sum((diff(diff(f[k](x), x), x))*sin(k*Pi*y), k = 1 .. infinity)+Sum(-f[k](x)*k^2*Pi^2*sin(k*Pi*y), k = 1 .. infinity))*sin(n*Pi*y), y = 0 .. 1) = Int((2*x^3+6*x*y*(1-y))*sin(n*Pi*y), y = 0 .. 1)

The infinite summations reduce to single terms as before.  The right-hand side may be evaluated explictly.  We this arrive at:

> 

ode := 1/2*(diff(f[k](x),x,x) - k^2*Pi^2*f[k](x)) = int(rhs(tmp)*sin(k*Pi*y), y=0..1) assuming k::integer;

(1/2)*(diff(diff(f[k](x), x), x))-(1/2)*k^2*Pi^2*f[k](x) = 2*(Pi^2*k^2*x^2+6)*x*(-(-1)^k+1)/(k^3*Pi^3)

Step 4:  We solve this ode along with the boundary conditions which we obtained in steps 1 and 2.

> 

dsol := dsolve({ode, ic1, ic2});

f[k](x) = -8*exp(k*Pi*x)*((-1)^k*Pi^2*k^2-k^2*Pi^2+6*(-1)^k-6)/(Pi^5*k^5*(exp(k*Pi)-exp(-k*Pi)))+8*exp(-k*Pi*x)*((-1)^k*Pi^2*k^2-k^2*Pi^2+6*(-1)^k-6)/(Pi^5*k^5*(exp(k*Pi)-exp(-k*Pi)))+4*((-1)^k-1)*x*(Pi^2*k^2*x^2+12)/(Pi^5*k^5)

Conclusion: The solution of the PDE is

> 

mysol := eval(sol, dsol);

u(x, y) = Sum((-8*exp(k*Pi*x)*((-1)^k*Pi^2*k^2-k^2*Pi^2+6*(-1)^k-6)/(Pi^5*k^5*(exp(k*Pi)-exp(-k*Pi)))+8*exp(-k*Pi*x)*((-1)^k*Pi^2*k^2-k^2*Pi^2+6*(-1)^k-6)/(Pi^5*k^5*(exp(k*Pi)-exp(-k*Pi)))+4*((-1)^k-1)*x*(Pi^2*k^2*x^2+12)/(Pi^5*k^5))*sin(k*Pi*y), k = 1 .. infinity)

> 

eval(rhs(mysol), infinity=30):
plot3d(%, x=0..1, y=0..1, style=patchcontour);


 

Download mw1.mw

Solution...

Rouben Rostamian 7711
May 07 2018
4 1
> 

restart;
> 

pde := diff(u(x,y),x,x) + diff(u(x,y),y,y)
          = 2*Pi*(2*Pi*y^2 - 2*Pi*y - 1) * exp(Pi*y*(1-y)) * sin(Pi*x);

diff(diff(u(x, y), x), x)+diff(diff(u(x, y), y), y) = 2*Pi*(2*Pi*y^2-2*Pi*y-1)*exp(Pi*y*(1-y))*sin(Pi*x)

> 

bc1 := u(0,y) = sin(Pi*y),  u(1,y) = exp(Pi)*sin(Pi*y);
bc2 := u(x,2) = exp(-2*Pi)*sin(Pi*x),  u(x,0) = u(x,1);

u(0, y) = sin(Pi*y), u(1, y) = exp(Pi)*sin(Pi*y)

u(x, 2) = exp(-2*Pi)*sin(Pi*x), u(x, 0) = u(x, 1)

The function "z(x.,y)=e^(Pi x) sin(Pi y)" satisfies the boundary conditions bc1, therefore

we change the independent variable from u to v through u(x, y) = z(x, y)+v(x, y),

Then v(x, y) would satisfy the homogeneous boundary conditions v(0, y) = 0, v(1, y) = 0.  

> 

z := (x,y) -> exp(Pi*x)*sin(Pi*y);

proc (x, y) options operator, arrow; exp(Pi*x)*sin(Pi*y) end proc

> 

pde_v := eval(pde, u(x,y) = z(x,y) + v(x,y));

diff(diff(v(x, y), x), x)+diff(diff(v(x, y), y), y) = 2*Pi*(2*Pi*y^2-2*Pi*y-1)*exp(Pi*y*(1-y))*sin(Pi*x)

In view of v's homogeneous boundary condition, we may express it in a Fourier series,
that is, an infinite sum of the form "(&sum;)`f__k`(y)*sin(k Pi x)."   However, the right-hand

side of the pde_v has only sin(Pi*x), therefore the infinite sum collapses to a single term.

Thus, we look for a solution of the pde_v in the form v(x, y) = f(y)*sin(Pi*x):

> 

eval(pde_v, v(x,y)=f(y)*sin(Pi*x)):
ODE := simplify(%/sin(Pi*x));

-f(y)*Pi^2+diff(diff(f(y), y), y) = 4*(-1/2+(y^2-y)*Pi)*Pi*exp(-Pi*y*(-1+y))

Thus the problem has been reduced to the solution of an ODE.

> 

dsol := dsolve(ODE);

f(y) = exp(Pi*y)*_C2+exp(-Pi*y)*_C1+exp(-Pi*y*(-1+y))

and we arrive at the solution of {pde + bc1}

> 

mysol := u(x,y) = eval(f(y), dsol) * sin(Pi*x) + z(x,y);

u(x, y) = (exp(Pi*y)*_C2+exp(-Pi*y)*_C1+exp(-Pi*y*(-1+y)))*sin(Pi*x)+exp(Pi*x)*sin(Pi*y)

We apply the boundary conditions bc2 to determine _C1 and _C2.  We have:

> 

simplify(eval(mysol, y=2) - bc2[1]);

0 = sin(Pi*x)*(exp(2*Pi)*_C2+exp(-2*Pi)*_C1)

For this to hold for all x, we need _C1 = _C2 = 0.  Thus mysol reduces to:

> 

mysol := eval(mysol, {_C1=0, _C2=0});

u(x, y) = exp(-Pi*y*(-1+y))*sin(Pi*x)+exp(Pi*x)*sin(Pi*y)

We have yet to apply the second condition in bc2, but there is no more freedom in

mysol.  But not to worry -- the problem is artificially concocted so that the second

condition in bc2 is automatically satisfied:

> 

eval(mysol, y=0) - eval(mysol, y=1);

u(x, 0)-u(x, 1) = 0

Conclusion:  The solution is:

> 

mysol;

u(x, y) = exp(-Pi*y*(-1+y))*sin(Pi*x)+exp(Pi*x)*sin(Pi*y)

> 

pdetest(mysol, {pde, bc1, bc2});

{0}

 

 

Download mw.mw

First 29 30 31 32 33 34 35 Last Page 31 of 53