tomleslie

Possibly something like...

Answered: tomleslie 13876

October 25 2018

1 4

the attached, although as VV has pointed you will have two watch for Maple's automatic simplification (and maybe other things I haven't considered!)

>

  restart;
  g:= proc(z)
           local f:=p->op([0..nops(p)],p):
           if   numelems(indets(z, name))=2
           then if   op(0,z)=`+`
                then if   member(-1, f~({f(z)}))
                     then return `-`;
                     else return `+`;
                     fi:
                elif op(0,z)=`*`
                then if   member(-1, f~({f(z)}))
                     then return `/`;
                     else return `*`;
                     fi:
                elif op(0,z)=`^`
                then return `^`;
                fi;
           else error "Does not have two operands":
           fi
      end proc:

  g(a+b);
  g(a-b);
  g(a*b);
  g(a/b);
  g(a^b);
  g(a*b+c);

Error, (in g) Does not have two operands

>

Download typrs.mw

Which is the 'variable'...

Answered: tomleslie 13876

October 24 2018

1 0

and which is a 'constant'. Unless you specify Maple (probably) just uses indets() which returns the indetermiinates in lexicographic order - which is alphabetic, except that all upper case names precede lower case ones.

Consider

indets(-v1+a);
indets(-V1+a);
{a, v1}
{V1, a}

so when you ask for

restart;
sign(-v1+a);
sign(-V1+a);

the first retuirns the sign of 'a', and the second returns the sign of V1. If you want too avoid this, specify the indeterminate, as in

sign(-v1+a, [v1]);
sign(-V1+a, [V1]);
-1
-1

which returns -1 in both cases

Something like?...

Answered: tomleslie 13876

October 24 2018

0 0

restart:
List1 := [v1,v2,v3,v4];
List2 := [a*b/c, a/c,a*b-2,d];
computation := subs( List1=~List2, v1*diff(expr1(x),x)+v3*diff(expr2(y),y)+ v3*v4);

Same graph - two ways...

Answered: tomleslie 13876

October 23 2018

1 3

In the attached your "pedal curve" is calculated two ways, with and without using a procedure.

Ouput "renders" better in a workeheet than it does on this site - honest

>

  restart;

  Trunc := proc(F::{`+`, procedure}, odr::posint := 2, v::(list(name)) := [x, y])
                 local f;
                 description " Truncates an algebraic equation to required degree";
                 f := `if`(F::procedure, F(v[]), F);
                 if not f::`+`
                 then error "Can't truncate 1-term expression";
                 else select(q->degree(q, v) <= odr, f);
                 end if;
           end proc:
  C := (x^2+y^2+12*x+9)^2-4*(2*x+3)^3:
  xp1 := 1:
  yp1 := 2:
  p21 := plots:-implicitplot(C,
                             x = -3 .. 3,
                             y = -5 .. 5,
                             colour = "Salmon",
                             gridrefine = 2,
                             size = [300, 300]
                            ):

  sol1 := [solve(C, y)]:
  f1 := expand
        ( eval
          ( C,
            [x = X+r, y = Y+s]
          )
        ):
  ltg1 := expand
          ( eval
            ( Trunc(f1, 1, [Y, X]),
              [X = x-r, Y = y-s]
            )
          ):
  lprp1 := -coeff(ltg1, y)*x+coeff(ltg1, x)*y+K1:
  for i to nops(sol1) do
      Pedal||i := rhs~( solve
                        ( eval
                          ( [ ltg1,
                              eval
                              ( lprp1,
                                [ K1 = solve
                                       ( eval
                                         ( lprp1,
                                           [x = xp1, y = yp1]
                                         ),
                                         K1
                                       )
                                ]
                              )
                            ],
                            s = eval(sol1[i], x = r)
                          ),
                          [x, y]
                        )[]
                      )[];
      p9||i := plot( [Pedal||i, r = -20 .. 20],
                     colour = "MediumSeaGreen"
                   );
  end do:
  plots:-display( p21, seq(p9 || i, i = 1 .. nops(sol1)),
                  size = [600, 600],
                  scaling = constrained,
                  caption = " Curve and it's Pedal curve"
                );

>

  restart;

  Trunc := proc(F::{`+`, procedure}, odr::posint := 2, v::(list(name)) := [x, y])
                 local f;
                 description " Truncates an algebraic equation to required degree";
                 f := `if`(F::procedure, F(v[]), F);
                 if not f::`+`
                 then error "Can't truncate 1-term expression";
                 else select(q->degree(q, v) <= odr, f);
                 end if;
           end proc:
  C := (x^2+y^2+12*x+9)^2-4*(2*x+3)^3:
  xp1 := 1:
  yp1 := 2:
  p21 := plots:-implicitplot(C,
                             x = -3 .. 3,
                             y = -5 .. 5,
                             colour = "Salmon",
                             gridrefine = 2,
                             size = [300, 300]
                            ):

  Pedal:=proc( Curv, xp, yp)
               local i, sol1, f1, ltg1,lprp1;
               sol1:= [ solve(Curv, y) ];
               f1:= expand
                    ( eval
                      ( Curv,
                        [x = X+r, y = Y+s]
                      )
                    );
               ltg1:= expand
                      ( eval
                        ( Trunc( f1, 1, [Y, X]),
                          [X = x-r, Y = y-s]
                        )
                      );
               lprp1 := -coeff(ltg1, y)*x+coeff(ltg1, x)*y+K1:
               for i to nops(sol1) do
                   Pedal||i := rhs~( solve
                                     ( eval
                                       ( [ ltg1,
                                           eval
                                           ( lprp1,
                                             [ K1 = solve
                                                    ( eval
                                                      ( lprp1,
                                                        [x = xp, y = yp]
                                                      ),
                                                      K1
                                                    )
                                            ]
                                           )
                                         ],
                                         s = eval(sol1[i], x = r)
                                       ),
                                       [x, y]
                                     )[]
                                   );
                    p9||i := plot( [Pedal||i[], r = -20 .. 20],
                                   colour = "MediumSeaGreen"
                   ):
               end do:
               return seq( p9||i, i=1..nops(sol1));
         end proc:

  Pedal2:=Pedal(C, xp1, yp1);

  plots:-display( [ p21, Pedal2 ],
                  size = [600, 600],
                  scaling = constrained,
                  caption = " Curve and it's Pedal curve"
                );

>

Download pedCur.mw

simplify(..,power)...

Answered: tomleslie 13876

October 22 2018

0 0

will work - most of the time. So for your two examples, see the attached

>	rand1:=rand(2..9): a:=rand1(); b:=rand1(); c:=rand1(); (asqrt(aa^b))^(1/c); simplify(%, power);

(1)

>	assume(x>0): sqrt(x^3*x^(3/4)); simplify(%,power);

(2)

>

Download pow.mw

Quick and dirty solution...

Answered: tomleslie 13876

October 21 2018

1 1

As has already been stsated the problem of fitting experimental data to a "model" is a whole subject on its own. A lot depends on how complicated you want to get! It is probably worth trying a simple approach first.

The attached file is concepyually about as simple as you can get, it simply determines the parameters which will minimize the 'penalty' function

sum( (x_model-x_exp)^2+(y_model-y_exp)^2)

where the sum is taken over all points. It would be fairly trivial to rewrite the 'fitter' procdure in the solution module for any other penalty function, or to add a second model function to the solution module.

The only real drawback whcih I had was that the Optimization:-NLPSolve() command, and I had to use the DirectSearch() package which is free to download from the MapleSoft Applications Centre, and *seems* gernally to be a bit more robust. This does however eman that you will be unable ro re-execute the attached unless you have the DirectSearch() package loaded.

>

restart;
#
# Define a solution module
#
  solmod:= module()
                  option package;
                  export getXPData, getSol, fitter:
                  local ModuleLoad, XPData, DEsol:
           #
           # Use the ModuleLoad() as an initalisation
           # to
           #
           # 1) read the experimental data from an
           #    external file. (OP will have to change
           #    the file path
           # 2) set up the solution procedure for the
           #    ODE system
           #
                  ModuleLoad:=proc()
                                   uses ExcelTools;
                                   local fname:="C:/Users/TomLeslie/Desktop/XPdata.xlsx",
                                         M:=Import(fname,1,"A2"),
                                   g:=9.81,
                                   m:=0.25,
                                   ode:= { D(x__1)(t) = x__2(t),
                                           D(x__2)(t) = -Drag*sqrt(x__2(t)^2+x__4(t)^2)*x__2(t)/m,
                                           D(x__3)(t) = x__4(t),
                                           D(x__4)(t) = -g-Drag*sqrt(x__2(t)^2+x__4(t)^2)*x__4(t)/m
                                         },
                                   ics:= { x__1(0) = 0,
                                           x__2(0) = v__0*cos(alpha),
                                           x__3(0) = 0,
                                           x__4(0) = v__0*sin(alpha)
                                         },
                                   solution:=dsolve(`union`(ode,ics),
                                                     numeric,
                                                     parameters=[alpha, v__0, Drag],
                                                     events=[[x__3(t), halt]]):
                                   M[1,1]:=0;
                                   return M, solution;
                            end proc:
           #
           # Execute the ModuleLoad() to return the experimental
           # data and the parameterized ODE system
           #

                XPData, DEsol:= ModuleLoad();
           #
           # Define a procedure which returns a matrix containing
           # t, x(t) and y(t), from the ODE system model, for any
           # supplied set of parameters alpha, velocity and Drag
           #
                getSol:=proc( p1, p2, p3 )
                              local tstep:=1/30, ans, data__halt:
                              interface(warnlevel=0):
                              DEsol(parameters=[p1*Pi/180,p2,p3]);
                              data__halt:=DEsol(10);
                              interface(warnlevel=3):
                              ans:=Matrix( [ seq( rhs~(DEsol(j)),
                                                  j=0..rhs(data__halt[1]), tstep
                                                )
                                           ]
                                         );
                              return Matrix( [ans[..,1], ans[..,2], ans[..,4]], scan=rows);
                       end proc;
           #
           # Define a procedure which simply returns the experimental
           # data matrix
           #
                getXPData:=proc()
                                return XPData;
                           end proc;
           #
           # Define a 'fitter' procedure which calculates a
           # simple least-square difference between the
           # experimental data and that returned by solving
           # the ODE system for any supplied set of parameters
           # alpha, velocity and Drag
           #
                fitter:=proc( p1, p2, p3 )
                              local M:=getSol(p1, p2, p3):

                              add( (M[j,2]-XPData[j,2])^2,
                                   j=1..min(op([1,1],M),op([1,1],XPData))
                                 )
                              +
                              add( (M[j,3]-XPData[j,3])^2,
                                   j=1..min(op([1,1],M),op([1,1],XPData))
                                 );
                        end proc;
         end module:

>

#
# Instantiate the module, load the plots() package and
# defimne a plot function which is going to be used a
# few times
#
  with(solmod):
  with(plots):
  doPlot:=(dat, col)->pointplot( dat[..,2], dat[..,3],
                                 style=point,
                                 symbol=solidcircle,
                                 symbolsize=16,
                                 labels=[x-position, y-position],
                                 labelfont=[times, bold, 20],
                                 color=col
                               ):

>

#
# Get the solution from the ODE system for the
# supplied set of paramters and generate a plot
# of y-position against x-position
#
  res:=getSol(47, 11.5, 0.8450638480e-2):
  pA:=doPlot(res, red):
#
# Get the experimental data and generate a plot
# of y-position against x-position
#
  XP:=getXPData():
  pB:=doPlot(XP, blue):
#
# Display the experimental and model data on the
# same graph
#
  display([pA,pB]);
#
# Close but no coconut!
#

>	# # Tried to use the 'fitter' procedure with NLPSolve # to find 'best' values for angle, velocity and # Drag, but this fails:-( # Optimization:-NLPSolve(fitter, 40..60, 10..15, 0.5e-02..1.0e-01);

Error, (in Optimization:-NLPSolve) no improved point could be found

>	# # Use the DirectSearch() package to obtain the best # fit. This throws a couple of warnings, but seems # to work # fitPars:=DirectSearch:-Search( fitter, [p1=40..60, p2=10..15, p3=0.5e-02..1.0e-01] );

Warning, initial point [p1 = .900000000000000, p2 = .900000000000000, p3 = .900000000000000] does not satisfy the inequality constraints; trying to find a feasible initial point

Warning, the new feasible initial point is [p1 = 46.7517737812633, p2 = 11.8882326721528, p3 = 0.979955380074218e-2]

[.19210558041243112, Vector(3, {(1) = 42.32561324590592, (2) = 11.574448061577609, (3) = 0.500000007987103e-2}), 207]

(1)

>

#
# Get the solution from the ODE with the values
# of parameters as returned by DirectSearch()
# and generate the plot of y-position versus
# x-position
#
  res:=getSol(convert(fitPars[2],list)[]):
  pA:=doPlot( res, red):
#
# Plot this solution agsainst that obtained for
# for the experimental data above
#
  plots:-display( [pA,pB],
                  title="Fitted Version",
                  titlefont=[times, bold, 20]);

>

Download projectile2.mw

The folly of 2D input...

Answered: tomleslie 13876

October 20 2018

1 1

If I check the "unknowns" prior to running the fsolve() command, then both 'f3' and 'Rsh' are reported as unknown. If I re-enter these with a simple copy+paste then hey presto, now they are 'known and the fsolve command works. Yet another clear cemonstration of twhy I wouldn't use 2D-input for a bet

If anyone can tell me the difference between the first and second assignments of 'Rsh' and 'f3' in the attached - I'd be vaguely interested. Or maybe I'll just stick to using 1D-input -which works

By the way, you don't really gain much for Digits settings higher than 20

>

restart; UseHardwareFloats := false; Digits := 40

(1)

>

f1 := Rs = -dvdi-Vt/(Io*exp(Voc/Vt))

>

f2 := Io = (Isc-Voc/Rsh)/(exp(Voc/Vt)-1)

>

fsolve([f1, f2, f3], fulldigits); `union`(`~`[indets]([f1, f2, f3], 'name')[])

Error, (in fsolve) number of equations, 3, does not match number of variables, 5

(2)

>

Rsh := 326; f3 := Impp = Isc-Io*exp((Impp*Rs+Vmpp)/Vt)-(Impp*Rs+Vmpp)/Rsh; `union`(`~`[indets]([f1, f2, f3], 'name')[]); fsolve([f1, f2, f3], fulldigits)

(3)

>

Download asoln.mw

The easy(?) way...

Answered: tomleslie 13876

October 19 2018

2 0

Your "model" is only valid until the projectile hits the ground - calculation beyond this point is meaningless.

You can use the 'events' option to stop the solution process when the y-position is zero.

It is then trivial to populate a matrix with all required data from t=0 to the 'halt' time, with any step size you want.

Use this 'model' data to compare with your 'experimental' data. Precisely how you do this comparison is a whole other subject!!

See the attached

Throw with air resistance

Constants

=

System of ODE

is position in x-direction

is the velocity in x-direction

is position in y-direction

is the velocity in y-direction

$odesystem := {ode1, ode2, ode3, ode4, x__1(0) = 0, x__2(0) = v__0*cos(alpha), x__3(0) = 0, x__4(0) = v__0*sin(alpha)}$

>

#
# Set up the ODE solution so that it "halts" when
# the projectile hits the ground - which occurs
# when the y-position ( ie x__3(t) ) is zero.
#
# The model is invalid beyond this time
#
  solution := dsolve(odesystem, numeric, events=[[x__3(t), halt]]):
  interface(warnlevel=0):
  data__halt:=solution(10);
  interface(warnlevel=3):
#
# Generate all data from this solution, in timesteps
# of tstep (OP seems to want 1/30) and store the
# results in a matrix whose columns are
#
# t x__1(t) x__2(t) x__3(t) x__4(t)
#
# Note that the time of hitting the ground may not be
# be a multiple of the chosen time step (1/30), so
# data for the halt time (probably) won't be included
# in the output matrix, but is available from the
# variable data__halt
#
# "Model" data from this matrix can then be compared with
# "experimental" data in a variety of ways
#
  tstep:=1/30:
  ans:=Matrix( [ seq( rhs~(solution(j)),
                      j=0..rhs(data__halt[1]), tstep
                    )
               ]
             );
#
# As an illustration, plot the vertical position against
# the horizontal position using the points from the above
# matrix
#
  plots:-pointplot( ans[..,2], ans[..,4],
                    style=point,
                    symbol=solidcircle,
                    symbolsize=16,
                    labels=[x-position, y-position],
                    labelfont=[times, bold, 20],
                    color=red
                  );
#
# In case it makes things easier to compare with 'experimental'
# data, extract clumns from the above complete matrix, to
# produce another matrix containing only time t, x-position
# x__1(t) and y-position x__3(t)
#
  ans2:= Matrix( [ans[..,1], ans[..,2], ans[..,4]], scan=rows);

data__halt := [t = 1.57327929242181, x__1(t) = 10.0194885388110, x__2(t) = 5.20670616099533, x__3(t) = -3.64291929955129*10^(-17), x__4(t) = -7.20221849558883]

(1)

>

Download projectile.mw

User protection...

Answered: tomleslie 13876

October 18 2018

0 0

MAple attempts to "protect" the user from the consequences of displaying *large* amounts of data. This "protection" can be overridden by setting the interface variable 'rtablesize', which is 10 by default. Although there is nothing wrong with mmcdara's solution, it will output any size of structured data - and if one of the dimensions happens to be 1000000, then you will get an awful lot of screenfuls of data. It is *probably* safer to stick with a "hard" limit as in the following where rtablesize is set to 30

>	restart; with( DEtools ): with( plots ): with( LinearAlgebra ): interface(rtablesize=30): ode2 := diff( x(t), t,t ) + diff( x(t), t ) - 2*x(t) = 0;

diff(diff(x(t), t), t)+diff(x(t), t)-2*x(t) = 0

(1)

>	ic2 := x(0)=X0, D(x)(0)=0; ic3 := eval(ic2,X0=3);

(2)

>	nsol := dsolve( {ode2,ic3}, x(t), type=numeric );

(3)

>

nsol(1);

(4)

>	eval(x(t), nsol(1));

(5)

>	odeplot( nsol, [[t,x(t)],[t,diff(x(t),t)]], 0..1, legend=[`x(t)`,`x'(t)`], title="Figure 1.1" );

>	odeplot( nsol, [t,x(t)-diff(x(t),t)], 0..10, title="Figure 1.2" );

>	Xeuler := dsolve( {ode2,ic3}, x(t), type=numeric, method=classical, stepsize=0.1 );

$proc (x_classical) 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_classical) else _xout := evalf(x_classical) end if; _dat := Array(1..4, {(1) = proc (_xin) local _xout, _n, _y0, _yn, _yl, _ctl, _octl, _reinit, _errcd, _fcn, _i, _yini, _pars, _ini, _par; option `Copyright (c) 2002 by the University of Waterloo. All rights reserved.`; table( [( "complex" ) = false ] ) _xout := _xin; _ctl := Array(1..9, {(1) = 2, (2) = 0., (3) = 0., (4) = .1, (5) = 50000, (6) = 0, (7) = 0, (8) = 1, (9) = 0}); _octl := Array(1..9, {(1) = 2, (2) = 0., (3) = 0., (4) = .1, (5) = 50000, (6) = 0, (7) = 0, (8) = 1, (9) = 0}); _n := trunc(_ctl[1]); _yini := Array(0..2, {(1) = 0., (2) = 3.}); _y0 := Array(0..2, {(1) = 0., (2) = 3.}); _yl := Array(1..2, {(1) = 0, (2) = 0}); _yn := Array(1..2, {(1) = 0, (2) = 0}); _fcn := proc (N, X, Y, YP) option `[Y[1] = x(t), Y[2] = diff(x(t),t)]`; YP[2] := -Y[2]+2*Y[1]; YP[1] := Y[2]; 0 end proc; _pars := []; if not type(_xout, 'numeric') then if member(_xout, ["start", "left", "right"]) then return _y0[0] elif _xout = "method" then return "classical" elif _xout = "numfun" then return round(_ctl[6]) elif _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "initial_and_parameters" then return [seq(_yini[_i], _i = 0 .. _n)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] elif _xout = "last" then if _ctl[3]-_y0[0] <> 0. then _xout := _ctl[3] elif _ctl[2]-_y0[0] <> 0. then _xout := _ctl[2] else error "no information is available on last computed point" end if elif _xout = "enginedata" then return eval(_octl, 1) elif _xout = "function" then return eval(_fcn, 1) 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, _yini) end if; if _ini <> [] then `dsolve/numeric/process_initial`(_n, _ini, _yini, _pars) end if; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if; _octl[2] := _y0[0]; _octl[3] := _y0[0]; for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; if _Env_smart_dsolve_numeric = true and type(_y0[0], 'numeric') then procname("right") := _y0[0]; procname("left") := _y0[0] end if; if _xout = "initial" then return [seq(_yini[_i], _i = 0 .. _n)] elif _xout = "parameters" then return [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] else return [seq(_yini[_i], _i = 0 .. _n), op(_pars)], [seq(_yini[_n+_i], _i = 1 .. nops(_pars))] end if else return "procname" end if end if; if _xout-_y0[0] = 0. then return [seq(_y0[_i], _i = 0 .. _n)] end if; _reinit := false; if _xin <> "last" then if 0 < 0 and `dsolve/numeric/checkglobals`(0, table( [ ] ), _pars, _n, _yini) then _reinit := true; if _pars <> [] then _par := {seq(rhs(_pars[_i]) = _yini[_n+_i], _i = 1 .. nops(_pars))}; for _i from 0 to _n do _y0[_i] := subs(_par, _yini[_i]) end do; for _i from _n+1 to _n+nops(_pars) do _y0[_i] := _yini[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do else for _i from 0 to _n do _y0[_i] := _yini[_i] end do end if end if; if _pars <> [] and select(type, {seq(_yini[_n+_i], _i = 1 .. nops(_pars))}, 'undefined') <> {} then error "parameters must be initialized before solution can be computed" end if end if; if not _reinit and _xout-_ctl[3] = 0. then [_ctl[3], seq(_yn[_i], _i = 1 .. _n)] elif not _reinit and _xout-_ctl[2] = 0. then [_ctl[2], seq(_yl[_i], _i = 1 .. _n)] else if _ctl[2]-_y0[0] = 0. or sign(_ctl[2]-_y0[0]) <> sign(_xout-_ctl[2]) or abs(_xout-_y0[0]) < abs(_xout-_ctl[2]) or _reinit then for _i to 9 do _ctl[_i] := _octl[_i] end do; for _i to _n+nops(_pars) do _yl[_i] := _y0[_i] end do; for _i to nops(_pars) do _yn[_n+_i] := _y0[_n+_i] end do end if; _ctl[3] := _xout; if Digits <= evalhf(Digits) then try _errcd := evalhf(`dsolve/numeric/classical`(_fcn, var(_ctl), _y0[0], var(_yl), var(_yn), Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), var(Array(1..3, 1..2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0, (3, 1) = 0, (3, 2) = 0})))) catch: userinfo(2, `dsolve/debug`, print(`Exception in classical:`, [lastexception])); if searchtext('evalhf', lastexception[2]) <> 0 or searchtext('real', lastexception[2]) <> 0 or searchtext('hardware', lastexception[2]) <> 0 then _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..3, 1..2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0, (3, 1) = 0, (3, 2) = 0})) else _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end if end try else try _errcd := `dsolve/numeric/classical`(_fcn, _ctl, _y0[0], _yl, _yn, Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..2, {(1) = 0, (2) = 0}), Array(1..3, 1..2, {(1, 1) = 0, (1, 2) = 0, (2, 1) = 0, (2, 2) = 0, (3, 1) = 0, (3, 2) = 0})) catch: _ctl[3] := _y0[0]; _errcd := [lastexception][2 .. -1] end try end if; if type(_errcd, 'list') or _errcd < 0 then userinfo(2, {dsolve, `dsolve/classical`}, `Last values returned:`); userinfo(2, {dsolve, `dsolve/classical`}, ` t =`, _ctl[2]); userinfo(2, {dsolve, `dsolve/classical`}, ` y =`, _yl[1]); for _i from 2 to _n do userinfo(2, {dsolve, `dsolve/classical`}, ` `, _yl[_i]) end do; if type(_errcd, 'list') then error op(_errcd) elif _errcd+1. = 0. then Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, maxfun limit exceeded (see ?dsolve,maxfun for details)", evalf[8](_ctl[3]) else Rounding := `if`(_y0[0] < _xout, -infinity, infinity); error "cannot evaluate the solution past %1, unknown error code returned from classical %2", evalf[8](_ctl[3]), trunc(_errcd) end if end if; if _Env_smart_dsolve_numeric = true then if _y0[0] < _xout and procname("right") < _xout then procname("right") := _xout elif _xout < _y0[0] and _xout < procname("left") then procname("left") := _xout end if end if; [_xout, seq(_yn[_i], _i = 1 .. _n)] end if end proc, (2) = Array(0..0, {}), (3) = [t, x(t), diff(x(t), 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_classical, ["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_classical, '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_classical, ["last", 'last', "initial", 'initial', "parameters", 'parameters', "initial_and_parameters", 'initial_and_parameters', NULL]) then _xout := convert(x_classical, '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_classical), 'string') = rhs(x_classical); 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_classical), 'string') = rhs(x_classical)) elif _xout = "solnprocedure" then return eval(_solnproc) elif _xout = "sysvars" then return _vars end if; if procname <> unknown then return ('procname')(x_classical) else _ndsol := 1; _ndsol := _ndsol; _ndsol := pointto(_dat[2][0]); return ('_ndsol')(x_classical) 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)

>	p1 := odeplot( nsol, [t,x(t)], 0..2, color=red ):

>	p2 := odeplot( Xeuler, [t,x(t)], 0..2, color=blue ):

>	display( [p1,p2], title="Figure 1.3" );

>	times := Array( [i/10$i=0..10, $2..10] );

Vector[row](20, {(1) = 0, (2) = 0, (3) = 0, (4) = 0, (5) = 0, (6) = 0, (7) = 0, (8) = 0, (9) = 0, (10) = 0, (11) = 0, (12) = 0, (13) = 0, (14) = 0, (15) = 0, (16) = 0, (17) = 0, (18) = 0, (19) = 0, (20) = 0})

(7)

>	t1 := dsolve( {ode2,ic3}, x(t), type=numeric, method=classical[foreuler], output=Array([i/10$i=0..10, $2..10]) );

(8)

>	t1[2,1][6,2];

(9)

>	t2 := dsolve( {ode2,ic3}, x(t), type=numeric, method=classical[abmoulton], output=Array([i/10$i=0..10, $2..10]) ):

>	header1 := <` t `\| `Forward Euler`\| `Predictor-Corrector`>:

>	header2 := <`---`\| `-------------`\| `-------------------`>:

>	< header1, header2, <Column(t1[2,1],1..2) \| Column(t2[2,1],2) > >;

(10)

>

Download opSize.mw

This seems to cover...

Answered: tomleslie 13876

October 17 2018

0 2

both of your examples. Not sure which one you actually want!?

restart;
sol1:=rsolve({U(n)=3*(U(n-3)+U(n-4))+6, U(0) = 6, U(1) = 6, U(2) = 6, U(3)=24}, {U}, 'makeproc'):
seq( sol1(j), j=0..10);
sol2:=rsolve({U(n)=U(n-1)+4*U(n-3) ,U(0) = 6,U(1) = 12,U(2) = 18, U(3)=42}, {U}, 'makeproc'):
seq( sol2(j), j=0..10);

Not completely clear (to me)...

Answered: tomleslie 13876

October 16 2018

0 0

precisely what you want, but maybe(?) something like

L:=[seq((L||j)(a_1, a_2,a_3), j=1..4)];
R:=[seq((R||j)(a_1, a_2,a_3), j=1..4)];
M:=Matrix(4,4, (i,j)->L[i]=R[j]);
M[2,3]:=eval( M[2,3], [a_1=1,a_2=3,a_3=4]):
M;

Simple addition...

Answered: tomleslie 13876

October 15 2018

0 0

my previous responses.

Although I note that you have (somehow?) deleted your original questions/worksheets. This just makes the current question almost impossible for any other responder to understand. Don't do this!

The attached will return the maximum entry in each column of the output table - although as I said in a previus response

stores the result in a table (which isn't necessarily the most convenient data structure for subsequent processing!)

>

  restart;
  with(Optimization):
  G:=a-b*p:
  z:=q-G:
  M:=int((z-u)*f(u), u=-3500..z):
  N:=int((u-z)*f(u), u=z..1500):
  beta:=alpha*K+(1-K)*r/p:
  E__pi:=(p-c)*(G+mu)-(c-v)*M-(p-c+s-beta*(p-r-c-d+s))*N:
  EXE__pi := eval( E__pi,
                   [ c = 35,
                     a = 100000,
                     b = 1500,
                     mu = -1000,
                     v = 10,
                     d = 3,
                     s = 3,
                     f(u)= exp((-1/800)*(u+1000)*(u+1000))/sqrt(2*3.14*800)
                   ]
                 ):
  DEP__pi:=diff(EXE__pi, p):
  DEQ__pi:=diff(EXE__pi, q):
  DER__pi:=diff(EXE__pi, r):
#
# Suppress printing of NLPSolve() output
#
# printf( "%2s %14s %13s %13s %13s\n\n",
#         "k", "maximum", "p", "q", "r"
#       );
#
# Set up output table
#
  myName:=table():
  for k from 0.1 by 0.02 to 0.40 do
      E1__pi:=eval(EXE__pi,[alpha = 0.5, K=k]):
      #printf("%a   ", k);
      try ans:= NLPSolve
                 ( E1__pi,
                   p = 40 .. 60,
                   q = 15000 .. 26000,
                   r = 0 .. 25,
                   initialpoint = {p = 50, q = 15000, r = 0},
                   maximize
                 ):
        #
        # Suppress printing of NLPSolve output
        #
        # printf( "%4.2f   %8.6e   %8.6e   %8.6e   %8.6e \n",
        #         k, ans[1], rhs~(ans[2])[]
        #      );
        #
        # Use outputs of NLPSolve to compute a
        # few things and store these in the
        # table myName, indexed by the value of
        # the loop variable k
        #
          myName[k] := eval( [ DEP__pi,
                               DEQ__pi,
                               DER__pi,
                               beta
                            ],
                            [ alpha = .5,
                              K = k,
                              ans[2][]
                            ]
                          ):
    catch "Error, (in Optimization:-NLPSolve) no improved point could be found":
        #
        # Print warning message when NLPSolve()
        # command fails
        #
          printf("\n**NLPSolve() failed when k=%4.2f**\n\n", k);
    end try:

  od:
#
# Check output table contents for a couple
# of index values
#
  myName[0.1];
  myName[0.28];
#
# Print the complete contents of the output
# table
#
  printf( "\n\n%2s%16s%17s%17s%17s\n\n",
          "k", 'DEP__pi', 'DEQ__pi', 'DER__pi', 'beta'
        );
  seq
  ( printf( "%4.2f%17.6e%17.6e%17.6e%17.6e \n",
            j, myName[j][]
          ),
    j in sort( [ indices(myName, 'nolist' ) ] )
  );

**NLPSolve() failed when k=0.20**

**NLPSolve() failed when k=0.36**

k         DEP__pi          DEQ__pi          DER__pi             beta

0.10    8.124235e-003    1.324435e-006    8.942607e-007    1.631055e-001
0.12    7.612606e-003    9.129944e-007   -1.005420e-006    1.650366e-001
0.14    7.942559e-002    4.729956e-005    1.622852e-004    1.669568e-001
0.16    3.069396e-001    1.945506e-004   -3.926488e-004    1.689244e-001
0.18   -4.559813e-002   -3.429199e-005   -2.476548e-005    1.708311e-001
0.22   -1.622543e-002   -1.491153e-005    6.926161e-005    1.746870e-001
0.24    1.663421e-003    5.050235e-008    4.267943e-006    1.766222e-001
0.26   -8.794443e-004   -4.779734e-006   -3.894070e-006    1.785537e-001
0.28   -1.131987e-001   -8.133517e-005    5.388903e-005    1.804810e-001
0.30    7.856244e-003    1.410988e-006   -6.168876e-005    1.824194e-001
0.32    6.452105e-001    4.175094e-004    1.243250e-004    1.843385e-001
0.34   -8.884629e-003   -6.830507e-006    1.669743e-005    1.862763e-001
0.38    2.951367e-002    1.540137e-005    1.012055e-005    1.901387e-001
0.40    5.357337e-003   -1.042016e-006   -1.274846e-001    2.000000e-001

>	# # Generate the maximum value in each column of the above # seq ( max ( seq ( myName[j][k], j in sort ( [ indices ( myName, 'nolist' ) ] ) ) ), k=1..4 );

(1)

>

Download optProb3.mw

Anything wrong with...

Answered: tomleslie 13876

October 15 2018

0 0

the attached?

Download vecDiff.mw

So many syntax errors...

Answered: tomleslie 13876

October 14 2018

0 2

that it is hard to know were to start!

Your code contains a mixture of "active" and "inert" commands: I am assumed that most of the latter should actually be the former
Multiple function definitions which do not depend on the passed argument - I mean why??
Use of the symbol for 'Pi' as an argument in many function definitions - just don't
A persistent pop-up for the interactive Plot Builder - why?

So I made a *lot* of guesses about what you actually intended, and came up with the attached. I have no doubt that some of my guesses may be wrong so you should read the attached very carefully.

The NLPsolve(..,maximise) command *failed* for a couple of values for the loop index. I didn't bother to work out why: just "filtered" these cases with a try-catch construct.

The achieved value of the maximum doesn't vary very much, and is always obtained for more or les the same values of the variables 'p' and 'q'. The value of the variable 'r' to produce the maximum does however vary with loop index.

>

restart;
with(Optimization):
G:=a-b*p:
z:=q-G:
M:=int((z-u)*f(u), u=-3500..z):
N:=int((u-z)*f(u), u=z..1500):
beta:=alpha*K+(1-K)*r/p:
E__pi:=(p-c)*(G+mu)-(c-v)*M-(p-c+s-beta*(p-r-c-d+s))*N:
EXE__pi := eval( E__pi,
                 [ c = 35,
                   a = 100000,
                   b = 1500,
                   mu = -1000,
                   v = 10,
                   d = 3,
                   s = 3,
                   f(u)= exp((-1/800)*(u+1000)*(u+1000))/sqrt(2*3.14*800)
                 ]
               ):
DEP__pi:=diff(EXE__pi, p):
DEQ__pi:=diff(EXE__pi, q):
DER__pi:=diff(EXE__pi, r):
printf( "%2s %14s %13s %13s %13s\n\n",
        "k", "maximum", "p", "q", "r"
      );
for k from 0.1 by 0.02 to 0.40 do
    E1__pi:=eval(EXE__pi,[alpha = 0.5, K=k]):
    #printf("%a   ", k);
    try ans:= NLPSolve
               ( E1__pi,
                 p = 40 .. 60,
                 q = 15000 .. 26000,
                 r = 0 .. 25,
                 initialpoint = {p = 50, q = 15000, r = 0},
                 maximize
               ):
         printf( "%4.2f   %8.6e   %8.6e   %8.6e   %8.6e \n",
                 k, ans[1], rhs~(ans[2])[]
               );
    catch "Error, (in Optimization:-NLPSolve) no improved point could be found":
          printf("\n**NLPSolve() failed when k=%4.2f**\n\n", k);
    end try:
    myName := table( eval( [ aa = DEP__pi,
                             bb = DEQ__pi,
                             cc = DER__pi,
                             dd = beta
                           ],
                           [ alpha = .5,
                             K = k,
                             p = 50.2439391403638,
                             q = 23142.7537807761,
                             r = 6.22630454447208
                           ]
                         )
                   ):
od:

k        maximum             p             q             r

0.10   3.601447e+05   5.049785e+001   2.324840e+004   6.346205e+000
0.12   3.601452e+05   5.049785e+001   2.324835e+004   6.027413e+000
0.14   3.601458e+05   5.049785e+001   2.324830e+004   5.693151e+000
0.16   3.601464e+05   5.049785e+001   2.324824e+004   5.345824e+000
0.18   3.601470e+05   5.049786e+001   2.324818e+004   4.977804e+000

**NLPSolve() failed when k=0.20**

0.22   3.601484e+05   5.049786e+001   2.324804e+004   4.187889e+000
0.24   3.601491e+05   5.049786e+001   2.324797e+004   3.762237e+000
0.26   3.601499e+05   5.049785e+001   2.324789e+004   3.313322e+000
0.28   3.601508e+05   5.049785e+001   2.324781e+004   2.839169e+000
0.30   3.601517e+05   5.049785e+001   2.324772e+004   2.338727e+000
0.32   3.601526e+05   5.049785e+001   2.324764e+004   1.807411e+000
0.34   3.601537e+05   5.049785e+001   2.324753e+004   1.245333e+000

**NLPSolve() failed when k=0.36**

0.38   3.601560e+05   5.049784e+001   2.324731e+004   1.129764e-002
0.40   3.601572e+05   5.049786e+001   2.324717e+004   0.000000e+000

>

Download optProb.mw

Well I don't understand the question...

Answered: tomleslie 13876

October 11 2018

0 8

but if I fix a few syntax errors, then I can plot X[i]/f(eta) for any value of X[i]. See the attached. Now what exactly is it that you want to know about this ode system?

As usual, plots look a lot more interesting in a Maple worksheet than they do when rendered on this site - II have no idea why this happens, but I encourage you to download and re-execute

>

  restart;
  with(plots):
  fcns := {T(eta), f(eta)}:
  ep:= .1: M := 1: kp := .5: n := 1:
  ec:= .1: pr := 1: s := .1: N := 5:
  sys:= diff(f(eta),eta$3)+f(eta)*diff(f(eta),eta$2)-diff(f(eta),eta)*diff(f(eta),eta)+ep*ep+(M+1/kp)*(ep-diff(f(eta),eta)) = 0,
        diff(T(eta),eta$2)+pr*(f(eta)*diff(T(eta),eta)-n*diff(f(eta),eta)*T(eta))+pr*(ec*diff(f(eta),eta$2)*diff(f(eta), eta$2)+ec*(M+1/kp)*diff(f(eta),eta)^2+s*T(eta)) = 0:
  bc:= f(0)=0, D(f)(0)=1, D(f)(N) = ep, T(0) = 1, T(N) = 0:
  R:= dsolve(eval({bc, sys}), numeric):
  psi := [-0.9e-1, -0.7e-1, -0.4e-1, -0.1e-1, 0, 0.1e-1, 0.4e-1, 0.7e-1, 0.9e-1]:
  for i to 9 do
      plt[i]:=odeplot( R, [eta, psi[i]/f(eta)], eta=0..5);
  end do:
  display( [seq(plt[i], i=1..9)]);

>

Download odeProb.mw

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

E-Mail Address:
Password:
Remember Me:	Automatically sign in on future visits

Ask a Question

Create a Post

13876 Reputation

20 Badges

MaplePrimes Activity

These are answers submitted by tomleslie

Possibly something like...

Which is the 'variable'...

Something like?...

Same graph - two ways...

simplify(..,power)...

Quick and dirty solution...

The folly of 2D input...

The easy(?) way...

User protection...

This seems to cover...

Not completely clear (to me)...

Simple addition...

Anything wrong with...

So many syntax errors...

Well I don't understand the question...

Save this setting as your default sorting preference?

Ask a Question

Create a Post

Generating PDF…

Save this setting as your default sorting preference?
Note: You can change your preference any time in your account settings
Don't show this again