acer

32405 Reputation

29 Badges

19 years, 350 days
Ontario, Canada

Social Networks and Content at Maplesoft.com

MaplePrimes Activity


These are answers submitted by acer

Your worksheet seems to have been last saved by you in Maple 13.

Using Maple 16.02 (the oldest I have at hand) it succeeds for your tau_l=0.8 if Digits is raised for the computations done prior to the problematic fsolve call, and if Digits is temporarily reduced for that call.

restart; with(plots)

Digits:=15:

The value of the parameters

omega_2 := .9892367:

eq_n := d-omega_4*omega_3/(omega_2*(omega_2*(1+n)^3+omega_2*(1+n)^2+omega_3*(1+n))):

fsolve({eq_n});

{n = .204400893732215}

x := 0.1e-1*T:

{nu = 1.03105988003254}

N1 := 1:

tau_c := 0.4551282e-1:

e1 := 9*(5.25*5)/(7*(24-10.35)*12):

tau_l := 0.8e-1:

eq_l3 := (1+n)*(1-tau_l)*(1-l2)/(Gamma*(1-tau_l-tau_upsilon)*(1-e1-l1))-(1-l3)/(nu*(1-l2));

1.22319829096998/Gamma-1.44893115172597+1.44893115172597*l3

{Gamma = 1.12560976598776, l3 = .250000000000000}

IOR := .2578598:

.83005166510992

.287475393180557

h := (Gamma*l1+l2)*h2/(1+n)+omega_3*h2*l3/(omega_2*nu*(1+n)^2);

.540810514922412

0.419629983047938e-1

.232264170334192

Upsilon := w*(nu*h2*(Gamma*l1+l2)*(1+n)*omega_2+omega_3*h2*l3)/(nu*((2+n)*(1+n)*omega_2+omega_3));

0.619481908316401e-1

cpp := (1/3)*theta_4*(2*(Upsilon-upsilon)+Gamma*w*h2*l1/(1+n)-upsilon);

0.112091197673597e-1

eq_chi := cpp*N4+(1+n)*omega_2*nu*exp(x)*chi = tau_upsilon*((Upsilon-upsilon)*omega_3/((1+n)^2*omega_2)+(Upsilon-upsilon)/(1+n)+Gamma*w*h2*l1/(1+n)-upsilon)+(1+r)*chi:

fsolve({eq_chi});

{chi = 0.530719545110145e-2}

oas := .522*cpp:

a := k-chi;

0.366558028536924e-1

EQ1 := a = a2/((1+n)*omega_2)+a3/(omega_2*(1+n)^2)+omega_3*a4/(omega_2^2*(1+n)^3):

The equations to solve

EQ2 := c = y-f*e1-(1+n)*omega_2*nu*exp(x)*k:

EQ4 := Gamma*(1-tau_l-tau_upsilon)*(1-e1-l1)/((1+n)*c1) = (1-tau_l)*(1-l2)/c2:

EQ6 := c2 = (1+(1-tau_a)*r)*c1/((1+rho)*nu*exp(x)):

EQ9 := (1+rho)*kappa_2*c2 = (1+n)*omega_2*c1:

EQ11 := Gamma*l1+(1-tau_l)*l2/(1-tau_l-tau_upsilon) = nu*(1+n)*(1/kappa_2-1)*(1+f*(1+n)/((1-tau_l-tau_upsilon)*w*Gamma*h2))/psi:

EQ12 := nu*exp(x)*a2 = (1-tau_l)*Gamma*w*h2*l1/(1+n)+epsilon_1+eta_1-tau_upsilon*(Gamma*w*h2*l1/(1+n)-upsilon)-(1+tau_c)*c1-f*e1:

EQ13 := nu*exp(x)*a3 = (1-tau_l)*w*h2*l2+(1+(1-tau_a)*r)*a2/omega_2+eta_2-tau_upsilon*(Upsilon-upsilon)-(1+tau_c)*c2-(1+n)*epsilon_1:

EQ14 := nu*exp(x)*a4 = (1-tau_l)*w*h2*l3/nu+(1+(1-tau_a)*r)*a3/omega_3+eta_3-tau_upsilon*(Upsilon-upsilon)-(1+tau_c)*c3:

EQ15 := (1+(1-tau_a)*r)*a4/omega_4+(1-tau_b)*(benef-gis)+gis = (1+tau_c)*c4:

EQ16 := tau_a*r*a+tau_c*c+tau_l*w*h+tau_b*(benef-gis)*N4 = eta_1*N1+eta_2*N2+eta_3*N3+(benef-cpp)*N4:

EQ17 := eta_2 = theta_2*(w*h2*l2+r*a2/omega_2):

EQ19 := f*e1/(c1+f*e1) = f_c1:

EQ21 := (tau_a*r*a+tau_l*w*h+tau_b*(benef-gis)*N4)/(r*a+w*h+(benef-gis)*N4) = tau:

EQs := {EQ10, EQ11, EQ12, EQ13, EQ14, EQ15, EQ16, EQ17, EQ18, EQ19, EQ2, EQ20, EQ21, EQ3, EQ4, EQ5, EQ6, EQ7, EQ8, EQ9}:

sv := {a2 = 0.1221680719e-1, a3 = 0.1221680719e-1, a4 = 0.1221680719e-1, benef = 0.1124557996e-1, c = .1412084225, c1 = 0.4267487699e-1, c2 = 0.5054102276e-1, c3 = 0.5985711412e-1, c4 = 0.6e-1, eta_2 = 0.4379459675e-2, eta_3 = 0.4e-2, f = 0.2810224635e-2, gis = 0.1e-2, kappa_2 = .7624704939, rho = .1236315697, tau_a = 0.1502430093e-1, tau_b = 0.7118933432e-1, theta_2 = 0.3946903560e-1, theta_3 = 0.39e-1, epsilon_1 = 0.2278732126e-1}:

Digits, oldDigits := 10, Digits;

10, 15

{a2 = 0.2181958727e-1, a3 = 0.1124241131e-1, a4 = 0.1842706633e-1, benef = 0.2325955123e-1, c = .1716741092, c1 = 0.5077912975e-1, c2 = 0.5518159721e-1, c3 = 0.5996575141e-1, c4 = 0.6516468394e-1, eta_2 = 0.1043313324e-1, eta_3 = 0.3091761489e-2, f = 0.3389901163e-2, gis = 0.6199270941e-2, kappa_2 = .8435924269, rho = .2996595542, tau_a = .1664750365, tau_b = -.6526622810, theta_2 = 0.9203906044e-1, theta_3 = 0.39e-1, epsilon_1 = 0.5153824973e-1}

15

eval(map(rhs-lhs, EQs), newsol);

{-0.153701e-9, -0.35502e-10, -0.21572e-10, -0.20677e-10, -0.148362e-10, -0.942622e-11, -0.93461e-11, -0.82466e-11, -0.30785e-11, -0.27822e-11, -0.14772e-11, 0.26576e-12, 0.105684e-11, 0.26276e-11, 0.33262e-11, 0.105704e-10, 0.114792e-10, 0.123763e-10, 0.137975e-10, 0.8559e-9}

NULL

Download ACCOLLEY_Delali_-_Demographics_3_ac.mw

restart;

 

`#mover(mo("="),mo("?"));`

`#mover(mo("="),mo("?"));`

 

`print/queryequal` := proc(a,b)
   uses Typesetting;
   mrow(mrow(Typeset(a),
        mo("⁢"),
        mover(mo("="),mo("?")),
        mo("⁢"),
        Typeset(b)));
end proc:

 

queryequal(X, Y);

queryequal(X, Y)

queryequal( sum(f(n),n=1..infinity), sqrt(Pi)/2);

queryequal(sum(f(n), n = 1 .. infinity), (1/2)*Pi^(1/2))

Download mover_ex.mw

In the last example above the result is an actual expression (which you could further manipulate programmatically if you wanted). It is an unevaluated function call to the unassigned global name queryequal.

In 2D Math entry mode you could get a similar rendering effect directly on input, by using the Layout palette for placing ? above =, etc. Here's what I got with that approach. I put in an extra space on either side of the symbol, and the result was this multiplication of terms. (I don't find that useful for further manipulation. But you could also just use the palette approach for non-executable 2D Math used for discourse rather than further computation.)

restart

sum((1/2)*f(n)*`#mover(mo("="),mo("?"))`*sqrt(Pi), n = 1 .. infinity)

sum((1/2)*f(n)*`#mover(mo("="),mo("?"))`*Pi^(1/2), n = 1 .. infinity)

lprint(%)

sum(1/2*f(n)*`#mover(mo("="),mo("?"))`*Pi^(1/2),n = 1 .. infinity)

Download mover_2dinput.mw

Your original expression assigned to r is or type `+`.

Hence when you map something over that the result will also be a call to `+`. The new summands will be whatever is returned by f, when passed the operands of r. That's some fundamental behaviour of map.

When your procedure f has a call to print as its last statement then it will return NULL, because print returns NULL. That's the difference. In your second example (using print) the result from your map call is NULL, as the addition of the various NULL return values from f.

When you have your procedure f return lists then they too will be added under `+`, when f gets mapped over r. But the manner in which you repeatedly concatenate global list L makes the lists be of different lengths, so adding them is not sensible.

Perhaps you would rather have the procedure f return the operand Z itself, or perhaps NULL explicitly forced, or perhaps Z^2 the result of the action.

restart;

f:=proc(Z) local res;
   global L;  
   res := Z^2;
   L:=[ op(L), res];
   print("L = ",L," Z=",Z);
   return res;  
end proc:

r:=2/(x^2+1)+1/(x^2+1)^2;

2/(x^2+1)+1/(x^2+1)^2

L:=[]:
map(Z->f(Z),r);

"L = ", [4/(x^2+1)^2], " Z=", 2/(x^2+1)

"L = ", [4/(x^2+1)^2, 1/(x^2+1)^4], " Z=", 1/(x^2+1)^2

4/(x^2+1)^2+1/(x^2+1)^4

L;

[4/(x^2+1)^2, 1/(x^2+1)^4]

different.mw

I leave it to you to decide what reasonable thing you want procedure f to return. (That changing list L is not such a choice...)

The 5rd bullet point in the Description of the Help page for map states,

   "The result of a call to map (unless the fold or reduce option
    described below is used) is a copy of expr with the ith
    operand of expr replaced by the result of applying fcn to the
    ith operand. This is done for all the operands of expr. For
    expressions of type atomic, map(fcn, expr) is identical to
    fcn(expr). For a table or rtable (Array, Matrix, or Vector), fcn
    is applied to each element instead of each operand."

So if you map f over a sum then the result is the sum of f applied to the operands (summands). If you map f over a product then the result is a product of f applied to the operands (multiplicands). And if you map f over another function call then the result is that same kind of call, now made over f applied to its original operands (arguments). Ie,

map(F, a*b*c/d);

  F(a)*F(b)*F(c)*F(1/d)

map(F, a+b+c);

      F(a)+F(b)+F(c)

map(F, g(a,b,c));

    g(F(a),F(b),F(c))

By the way, I don't see why you want to build list L in such a recursive (augmenting) manner. More efficient could be something like, say,

r:=2/(x^2+1)+1/(x^2+1)^2;

2/(x^2+1)+1/(x^2+1)^2

[seq(Z^2, Z=[op(r)])];

[4/(x^2+1)^2, 1/(x^2+1)^4]

map(Z->Z^2, [op(r)]);

[4/(x^2+1)^2, 1/(x^2+1)^4]

I prefer the first approach below (as Tom suggested) since it's a very easy way to avoid the inefficiency of recomputing h(a) multiple times.

restart;

ode := diff(f(x),x) = f(x):

initial_condition := f(0) = 1:

a := 1:

 

# As Tom suggested

hlp := dsolve({ode, initial_condition}, numeric,
              output = listprocedure):
F := eval(f(x), hlp):

 

fsolve(F(x) - F(a) = 1);

1.313261761

# Alternatively

Fa := F(a):
fsolve(x -> F(x) - Fa - 1);

1.313261761

# And without listprocedure

h := dsolve({ode, initial_condition}, numeric):

 

ha := rhs(h(a)[2]):
fsolve(x -> rhs(h(x)[2]) - ha - 1);

1.313261761

ha := eval(f(:-x),h(a)):
fsolve(x -> eval(f(:-x),h(x)) - ha - 1);

1.313261761

Download ds_rf.mw

Yes, a 3D space-curve can be colored.

(I dug this out of an old directory in which I had worksheets that colored space-curves of knots, eg. by torsion and curvature. The procedure SC constructs a space-curve substructure as well as its corresponding color substructure. Apart from adjusting the cfunc for your example I had to make no other changes. That's because coloring a space-curve is a pretty basic goal. I have previously submitted a request for it in the plots:-spacecurve command.)

I prefer this over using tubeplot -- at least for my purposes -- for aesthetic, performance, and memory considerations.

restart:

kernelopts(version);

`Maple 2015.2, X86 64 LINUX, Dec 20 2015, Build ID 1097895`

sys := {diff(x(t),t)=v(t), diff(v(t), t)=cos(t), x(0)=-1, v(0)=0}:

sol := dsolve(sys, numeric, output=listprocedure):
X,V := eval([x(t),v(t)], sol)[]:

 

cfunc:=proc(tt) local xt,vt;
  xt,vt := X(tt),V(tt);
  ColorTools:-Color(piecewise(xt>=0 and vt>=0, "Green",
                              xt>=0, "Red",
                              vt>=0, "Blue", "Gold"))[]:
end proc:

 

SC:=proc( T::list, trng::name=range(realcons),
          {numpoints::posint:=3, colorfunc::procedure:=NULL} )

  local M, C, t, a, b, i;

  t := lhs(trng);

  a,b := op(evalf(rhs(trng)));

  M := [seq(eval(T,t=a+i*(b-a)/(numpoints-1)),i=0..numpoints-1)];

  if colorfunc=NULL then C := NULL; else

    C := :-COLOUR(':-RGB',hfarray([seq(colorfunc(a+(i-1)*(b-a)/(numpoints-1)),
                                       i=1..numpoints)]));

  end if;

  plots:-display(:-PLOT3D(:-CURVES(M,C),':-THICKNESS'(5)), _rest);

end proc:

 

CodeTools:-Usage(
  SC( [t,X(t),V(t)], t=0..4*Pi, numpoints=200, colorfunc=cfunc,
      labels=[t,x(t),v(t)]) );

memory used=8.07MiB, alloc change=32.00MiB, cpu time=101.00ms, real time=102.00ms, gc time=0ns

 

Download spacecurve_shaded.mw

I used an hfarray for some backwards compatibility, but in more recent Maple I could have used a Array with datatype=float[8].

I didn't use the px(t) and pv(t) and their piecewises, only because they weren't strictly necessary. I'm sure that you could see how they might be incorporated instead.

note: If you don't need to alter Digits while calling these, you can squeeze out faster performance by adding these lines after extracting X and V from sol.
   X:=subsop(3=[remember,system][],eval(X)):
   V:=subsop(3=[remember,system][],eval(V)):

This memoization improves performance in such code as this which invokes X and V more than once at the same t values.

edit: I'll mention that I deliberately made cfunc less efficient than necessary. It could instead have the RGB values (range 0..1) hard coded in the piecewise. Instead it needlessly calls Color and rips those values out, each time it is called. I just wanted to illustrate using the colors by name. So it could be faster, though for this example I expect the DE computation to dominate the cost of color generation. Let me know if you'd like me to show such an edit. The particular RGB values -- or the piecewise -- could even be substituted programmatically into the cfunc proc (once, up front). You could also call that cfunc yourself, and see what it returns directly.

A different and efficient approach for SC would be to have it call `spacecurve` or odeplot itself, rip out the numeric data from its resulting CURVES substructure, and then compute the COLOR substructure using that numeric data directly. That is how Library improved revisions to those routines would ideally be extended to implement the usual coloring options.

The angle-bracket shortcut syntax for Matrix construction supports this functionality. For example,

a := Matrix(2,3,symbol=A);

a := Matrix(2, 3, {(1, 1) = A[1, 1], (1, 2) = A[1, 2], (1, 3) = A[1, 3], (2, 1) = A[2, 1], (2, 2) = A[2, 2], (2, 3) = A[2, 3]})

b := Matrix(3,3,symbol=B);

b := Matrix(3, 3, {(1, 1) = B[1, 1], (1, 2) = B[1, 2], (1, 3) = B[1, 3], (2, 1) = B[2, 1], (2, 2) = B[2, 2], (2, 3) = B[2, 3], (3, 1) = B[3, 1], (3, 2) = B[3, 2], (3, 3) = B[3, 3]})

c := Matrix(2,2,symbol=C);

c := Matrix(2, 2, {(1, 1) = C[1, 1], (1, 2) = C[1, 2], (2, 1) = C[2, 1], (2, 2) = C[2, 2]})

<a|c>;

Matrix([[A[1, 1], A[1, 2], A[1, 3], C[1, 1], C[1, 2]], [A[2, 1], A[2, 2], A[2, 3], C[2, 1], C[2, 2]]])

<a,b>;

Matrix([[A[1, 1], A[1, 2], A[1, 3]], [A[2, 1], A[2, 2], A[2, 3]], [B[1, 1], B[1, 2], B[1, 3]], [B[2, 1], B[2, 2], B[2, 3]], [B[3, 1], B[3, 2], B[3, 3]]])

Download Matrix_stack_augment.mw

That syntax is so simple that having dedicated commands for it would be overkill.

You may also use this to stack/augment a Matrix by a Vector (of the appropriate orientation).

You can suppress the visual display of results by using a fill colon (instead of semicolon) as the statement terminator.

For sorting there are various kinds of example. You should explain the specific goals, explicitly.

You can pass the extra option  output=string  to the CodeGeneration:-Matlab command.

Eg,
   CodeGeneration:-Matlab(func, output=string);

Then you can write that string to a file, using fprintf, or FileTools:-Text:-WriteString, etc.

It sounds like you are looking for a 2D parametric plot.

For example (using made up values for the other parameters),

plot([eval([x(t),y(t)],[A=0.9,v=2,m=1,n=1,alpha=1,beta=1])[], t=0.3..1.5],
        labels=["x","y"]);

With the code you supplied, the following both work:

  ribbonplot5([cos, sin, cos + sin], -Pi .. Pi,numpoints=20);
  ribbonplot5([cos(x), sin(x), cos(x) + sin(x)],x=-Pi .. Pi,numpoints=20);

Consider this simpler plotting example, and these two different kinds of calling sequence,

  plot(sin, 0..Pi);
  plot(sin(x), x=0..Pi);

The first uses a so-called operator-form calling calling sequence syntax, and the second uses a so-called expression-form calling sequence. (Also see the calling sequences shown on this Help page.)

The first passes an operator (procedure), which the second passes an expression with an unknown (here, the name x).

It does not make good sense to mix and match those two distinct kinds of usage. It does not make good sense to pass an operator (procedure) as the first argument and a name in the second argument. That would be muddled up and poor programming.

That ribbonplot5 procedure simply doesn't do as good a job as plot at checking for such a mismatch and producing a helpful error message. Eg,

  plot([sin, cos], x=0..Pi);
Error, (in plot) expected a range but received x = 0 .. Pi

note: This kind of distinction between expression-form and operator-form calling sequence is present elsewhere, eg. for fsolve, evalf(Int(...)), Optimization:-Minimize, etc.

You might use the orientation 3D plotting option.

For example,

restart:

with(plots): with(plottools):

d := 1.25:
a := 1:

chl := [[a/2, a/2, a/2],
          [a/2, -a/2, -a/2],
          [-a/2, a/2, -a/2],
          [-a/2, -a/2, a/2]]:
f := [[[-d/2, -d/2, -d/2], [-d/2, d/2, -d/2], [d/2, d/2, -d/2], [d/2, -d/2, -d/2], [-d/2, -d/2, -d/2]],
      [[-d/2, -d/2, d/2], [-d/2, d/2, d/2], [d/2, d/2, d/2], [d/2, -d/2, d/2], [-d/2, -d/2, d/2]],
      [[-d/2, d/2, d/2], [-d/2, d/2, -d/2], [d/2, d/2, -d/2], [d/2, d/2, d/2], [-d/2, d/2, d/2]],
      [[-d/2, -d/2, d/2], [-d/2, -d/2, -d/2], [d/2, -d/2, -d/2], [d/2, -d/2, d/2], [-d/2, -d/2, d/2]]]:
l := seq(seq(line(op(i-1,j),
              op(i, j),
                        colour = black,
                 style = surface),
            i= 2..5), j in f):    
t := seq(line([0,0,0],
            i,
            colour = black,
            thickness = 30,
            style = surface),
            i in chl):
hydr := seq(sphere(3*i/4,
                   a/4,
                   colour = gray,
                   style = surface),
                   i in chl):
char := sphere([0,0,0],
               a/3,
               colour=black,
               style = surface):

display(l, t, hydr, char, axes = none, orientation=[45,90,0]);

display(l, t, hydr, char, axes = none, orientation=[45,54,0]);

display(l, t, hydr, char, axes = none, orientation=[135,90,0]);

Download tareaclase6_ac.mw

An imported image (as float[8] Array) can be used for the coloring in a density-plot structure.

Once appropriately scaled, this can be translated to arbitrary x-y positions.

The following code could be made more efficient (and easy, with a re-usable procedure to scale/translate/etc). But I'd like to know whether it will suffice, before making the effort.

Naturally one can change the Read line below, to instead use some other external image file. (The OP didn't provide an image file.) All the code below needs is the usefully scaled (70x70x3 works ok) Array assigned to img.

restart;

with(plottools, transform): with(ImageTools):

img := Scale(Read(cat(currentdir(),"/smdrone.jpg")),1..75)[1..70]:

R := Rotate(img,':-right'):

nc, nr := LinearAlgebra:-Dimensions(R[..,..,1]);

70, 70

Q:=PLOT(GRID(-0.15..0.15,-0.15..0.15,
             Matrix(nc,nr,'datatype'=':-float[8]'),
             COLOR(RGB,R)),STYLE(PATCHNOGRID)):

plots:-display(
  transform((x,y)->[x+Pi/2,y+1.1])(Q),
  transform((x,y)->[x+Pi/5,y+0.7])(Q),
  transform((x,y)->[x+5*Pi/6,y+0.7])(Q),
  transform((x,y)->[x+Pi+0.2,y+0.1])(Q),
  transform((x,y)->[x-0.2,y+0.1])(Q),

  plot(sin(x),x=0..Pi),

  scaling=constrained, view=[-0.4..Pi+0.4,0..1.5],
  axis[1]=[tickmarks=[seq(i*Pi/4=i*Pi/4,i=0..4)]],
  size=[500,350]);

Download image_plot.mw

I suspect that the GUI will be generally less sluggish using such density-plots than it would with, say, colored point-plot assemblies.

That density-plot structure could actually be formed via the plots:-densityplot command, using the Array from the image. I simply found it easier to form it directly, as above.

The Fractals:-EscapeTime commands produce images (in the float[8] rtable sense), not plots.

You can "animate" a discrete number of images (or even calls to the Mandelbrot command at modest size, it's that fast) using Explore.

But you cannot nicely export that animated rendering to an external file. Maple's directly export functionality is limited to exporting a sequence of actual plots. Converting those images (rtables) into, say, densityplots would produce results too large and unwieldy for the GUI and its export drivers to handle well (if at all).

What you could do, instead, is export the images (rtables) directly to separate external image files (.gif, .png, what have you), using the ImageTools:-Write command. And then use an external 3rd party tool on your OS to assemble a movie format file (mpeg2, mp4, .gif, etc) from those.

Here are two alternatives to plots:-display, one using DocumentTools:-Tabulate and the other using DocumentTools primitives.

If you really want you could use more Table nesting (and suppressed borders) to get precisely the layout you showed. Perhaps more effort than it's worth...

(Sorry, I don't think that this forum will properly inline this worksheet here.)

question5_ac.mw

I don't why you wouldn't expect a(i)=0 and c(i)=0 to follow from your given a(0)=0 and c(0)=0. Perhaps you (or I) have made a mistake?

For fun, procedures Ap and Cp below take a0,c0 initial values as 2nd and 3rd arguments respectively.

The procedures Ap and Cp are programmatically generated from the recurrence equations. (I didn't feel like wrestling with rsolve.) So you can regenerate them for some modest edits of the equations.

In the example I use a0=0.1 and c0=0.1.

restart

eqns := {a(i) = (1+10^(-5)*(0.1e-1*(i-1)-0.7e-2)/(0.1e-3))*a(i-1)+.1*10^(-5)*c(i-1), c(i) = 0.1e-1*10^(-5)*(i-1)*a(i-1)/(0.1e-3)+((1-10^(-5))*.1)*c(i-1)}

{a(i) = (.9983000000+0.1000000000e-2*i)*a(i-1)+0.1000000000e-5*c(i-1), c(i) = 0.1000000000e-2*(i-1)*a(i-1)+0.9999900000e-1*c(i-1)}

Ap := subs(_d = subs([a = Ap, c = Cp], subsindets(eval(a(i), eqns), specfunc({a, c}), proc (u) options operator, arrow; (op(0, u))(op(u), a0, c0) end proc)), proc (i, a0, c0) option remember; if not i::nonnegint then return ('procname')(i, a0, c0) elif i = 0 then a0 else _d end if end proc)

proc (i, a0, c0) option remember; if not i::nonnegint then return ('procname')(i, a0, c0) elif i = 0 then a0 else (.9983000000+0.1000000000e-2*i)*Ap(i-1, a0, c0)+0.1000000000e-5*Cp(i-1, a0, c0) end if end proc

Cp := subs(_d = subs([a = Ap, c = Cp], subsindets(eval(c(i), eqns), specfunc({a, c}), proc (u) options operator, arrow; (op(0, u))(op(u), a0, c0) end proc)), proc (i, a0, c0) option remember; if not i::nonnegint then return ('procname')(i, a0, c0) elif i = 0 then c0 else _d end if end proc)

proc (i, a0, c0) option remember; if not i::nonnegint then return ('procname')(i, a0, c0) elif i = 0 then c0 else 0.1000000000e-2*(i-1)*Ap(i-1, a0, c0)+0.9999900000e-1*Cp(i-1, a0, c0) end if end proc

``

plots:-display(plots:-listplot([seq(Ap(i, .1, .1), i = 1 .. 50)], color = red, style = point), plots:-listplot([seq(Cp(i, .1, .1), i = 1 .. 50)], color = blue, style = point))

``

Download rsolve_for_system_of_recursive_equations_acc.mw

There are various ways to programmatically make the result appear in terms of J/mol as units. See attachment below for three of them.

It numerically better to round final results to 3 digits instead of calling evalf(...,3) or evalf[3](...) which do computations at lower precision and can incur unnecessarily greater roundoff error than usual.

I am deliberately avoiding using right-click units-formatting or numeric (float) formatting via the GUI.

(For some reason this forum is not letting me inline my revised Document, sorry).
Download question4_ac.mw

In terms of programmatically determining significant digits you might try either the Tolerances or the ScientificErrorAnalysis packages, depending on your concept of error in the supplied input float values.

First 75 76 77 78 79 80 81 Last Page 77 of 336