tomleslie

It is pretty...

Answered: tomleslie 13876

November 25 2022

0 7

trivial to plot the analytical solution for the same parameters as the numerical solutions - either separately or together, see the attached.

One would have to say that the agreement is terrible - a typo somewhere maybe?

>

`λ__vals` := [0.96e-2, 0.10e-2, 0.93e-2, 0.75e-2, 0.31e-2]; k__vals := [.1, .3, .5, .7, .9]

>

colors := [black, red, blue, green, orange, cyan, purple, yellow, navy]; for j to numelems(k__vals) do sol := dsolve(eval([ode, bc], [k = k__vals[j], We = .1, Br = .3, x = 0, Q = Q__vals[j], lambda = `λ__vals`[j]]), numeric); pltn[j] := plots:-odeplot(sol, [y, theta(y)], y = 0 .. 1, axes = boxed, color = colors[j]) end do; plots:-display([seq(pltn[j], j = 1 .. 5)], title = "Numerical Solution")

>

Analytical_sol := -(3*y+3+3*x^2*(1/2))*((-2/3+(-(5/12)*k^2+(1/6)*k-5/12)*Br)*(1+(1/2)*x^2)^5-lambda*Br*((-(1/12)*k^2+(1/6)*k-3/4)*y+Q*(k-1))*(1+(1/2)*x^2)^4-(((1/6)*k+7/6)*y+Q)*Br*((-(1/2)*k+1/2)*y+Q)*(1+(1/2)*x^2)^3+We*Br*((1/4)*(k-1)^2*y^2-(1/3)*Q*(k-5)*y+Q^2)*y*(1+(1/2)*x^2)^2-((1-k)*y+Q)*Br*y^2*Q*(1+(1/2)*x^2)+Q^2*y^3*Br)/(4*(1+(1/2)*x^2)^6)

>

Download plots.mw

It would be helpful...

Answered: tomleslie 13876

November 24 2022

0 0

if you posted the complete worksheet for your calculation using the big green up-arrow in the Mapleprimes toolbar.

In the attached I have performed the calculation you seem to want in two completely different ways.

The first uses Maple's bult-in numapprox() and orthopoly() packages. The second uses an "add-on" package (which you probably don't have?) called OrthogonalExpansions() which is available from the Maple Application Centre. These return the same answers: the desired series is

-0.577350269189626*T(0, x) - 0.309401076758503*T(1, x) - 0.0829037686547607*T(2, x)

The coefficients in this series disagree with those which you report - but since you do not specify how you obtained your coefficients, there isn't much I can say?

When converted to a simple polynomial in 'x', both give

-0.494446500534865 - 0.309401076758503*x - 0.165807537309521*x^2

and both return your measure of "error" as 0.02776926226

>

  restart;
#
# Use built-in approach with numapproz and orthopoly packages
#
# NB numapprox(0 will only give floating point coefficients
# for the series. Retyrned series is claimed to be accurate
# to 10^-(Digits) - Diigits is on default, ie 10
#
  u:=1/(x-2);
  ser0:= numapprox:-chebyshev(u, x=-1..1):
#
# Restrict to first three terms, and use orthopoly() to
# evaluate the ChebyshevT polynomials. Note that the latter
# returns the same as using the OrthogonalExpansions package
# below
#
  ser2:=add(op([1..3], ser0));
  ser3:=eval( ser2, T=orthopoly[T]);
#
# Compute OP's definition of "error"
#
  evalf(Int(abs(u - ser3), x = -1 .. 1));

(1)

>

  restart;
  with(OrthogonalExpansions):
  u:=1/(x-2);
#
# Generate ChebyshevT series of 'u' above over the
# range x=-1..1
#
  ser := ChebyshevTSeries(u, x = -1 .. 1, 'n', 'Coefficients');
#
# Restrict series to th first three terms, ie, ChebyshevT(0, x)
# ChebyshevT(1, x), and ChebyshevT(2, x)
#
# Return the restricted series and the coefficients. Compare the
# latter with the values in 'ser2' above
#
  n:= 2;
  ser2:=value(ser);
  evalf(value(Coefficients));
#
# Use the orthopoly() package to convert the ChebyshevT polynomials
# Compare 'ser3' below with 'ser3' in previous execution group
#
  ser3:=evalf(eval(ser2, [ seq( ChebyshevT(i, x)=orthopoly:-T(i,x), i=0..2)]));
#
# Check OP's definition of the "error" between the series
# expansion and the orifginal function
#
  evalf(Int(abs(u - ser3), x = -1 .. 1));

Sum(piecewise(i = 0, 1/Pi, 2/Pi)*(int(ChebyshevT(i, x)/((x-2)*(-x^2+1)^(1/2)), x = -1 .. 1))*ChebyshevT(i, x), i = 0 .. n)

-(1/3)*3^(1/2)*ChebyshevT(0, x)+2*(Pi-(2/3)*Pi*3^(1/2))*ChebyshevT(1, x)/Pi+2*(4*Pi-(7/3)*Pi*3^(1/2))*ChebyshevT(2, x)/Pi

(2)

>

Download cheby.mw

Consider the "toy" example attached...

Answered: tomleslie 13876

November 22 2022

1 2

This illustrates how I'd probably do it.

The ODE system is a simple predator-prey model.
A list of solution matrices is obtained by varying one of the initial conditions
The list of solution matrices is plotted as an animation by using the display() command with the option, insequence=true
The list of solution matrices is exported to a (matlab) binary file - I think this is about the only option if one wants to export a list of matrices (as distinct from a single matrix). Obviously you will have to change the path name for this file to something appropriate for your machine
restart;
Import the list of solution matrices generated by 4 above
Repeat 3 above

>

  restart;
  with(plots):
  odesys:= [ diff(x(t),t) = x(t) - x(t)*y(t),
             diff(y(t),t) = x(t)*y(t) - y(t)
           ] :
  tv:= Array([seq( j, j=0..10,0.1)]):
  ics:= [x(0)=k, y(0)=1]:
#
# Solve odesys, for different initial conditions
# and extract the matrix which contains the results
# to be plotted
#
# 'sol' is a list of such matrices
#
  sol:= [ seq( dsolve( [odesys[], eval(ics, k=j)[]],
                       numeric,
                       output=tv)[2,1][..,2..3],
               j=1.5..2.5, 0.1
             )
        ]:
#
# Display the results as an animation
#
  display( [seq( listplot( sol[j]), j=1..numelems(sol))], insequence=true);
#
# Export the list of solution matrices
#
  ExportMatrix( "C:/Users/TomLeslie/Desktop/test.mat", sol, target=Matlab);

(1)

>	restart; with(plots): # # Import the list of solution matrices which was generated above # ans:=[ImportMatrix( "C:/Users/TomLeslie/Desktop/test.mat", output=matrices)]: # # Display the results as an animation # display( [seq( listplot( ans[j]), j=1..numelems(ans))], insequence=true);

>

Download toyAnim.mw

Ah - the Lotka-Volterra predator-prey sy...

Answered: tomleslie 13876

November 21 2022

0 0

This is actually used as an example (with multiple plotting possibilities) on the help page for DEPlot() .

There are *a lot* of options for this command. One possibility for your specific set of equations is shown in the attached.

If you don't want actual solution curves (for which I admit I just made up initial values), then you can just delete the lines in the attached

[ [x(0)=2,   y(0)=1],
  [x(0)=0.5, y(0)=3],
  [x(0)=1.5, y(0)=1.5]
],
linecolor=[ green, orange, yellow],

>

  DETools:-DEplot( [ diff(x(t),t) = x(t) - x(t)*y(t),
                     diff(y(t),t) = x(t)*y(t) - y(t)
                   ],
                   [x(t),y(t)],
                   t=0..10,
                   x=0..3,
                   y=0..3,
                   [ [x(0)=2,   y(0)=1],
                     [x(0)=0.5, y(0)=3],
                     [x(0)=1.5, y(0)=1.5]
                   ],
                   linecolor=[ green, orange, yellow],
                   arrows=mediumfill,
                   color=magnitude,
                   dirfield=[10,10],
                   axes=boxed,
                   labels=["prey/density", "predator/density"],
                   labelfont=[times, bold, 12]
                 );

>

Download predprey.mw

It depends...

Answered: tomleslie 13876

November 21 2022

0 0

on whether you want a 2-D contour plot or a 3-D contour plot! Both come with many options. I suggest you read the help at ?contoutplot.

Examples of 2D and 3D plots for your function are shown in the attached

>

  restart;
  with(plots):
  f:=(x,y)-> (x - 3)*(y - 2) - 2;
  contourplot( f(x,y),
               x=-10..10,
               y=-10..10,
               contours=20,
               filledregions=true,
               coloring=[red, blue]
             );
  contourplot3d( f(x,y),
                 x=-10..10,
                 y=-10..10,
                 contours=20,
                 filledregions=true,
                 coloring=[red, blue]
               );

>

Download conPlot.mw

Two points...

Answered: tomleslie 13876

November 20 2022

0 0

Error is not new - it occurs in Maple 2021.2. (i haven't bothered to check anything earlier
"too many levels of recursion" is untrappable. see the excerpt from the Maple help below, with emphasis added

Some exceptions are untrappable. This means that a catch clause will never match such an exception, even if the catchString is missing. The finalStatSeq is also not executed if an untrappable exceptions occurs. The following exceptions are untrappable:
• Computation timed out (this can only be caught by timelimit, which raises a "time expired" exception, which can be caught).
• Computation interrupted (user pressed Ctrl+C, Break, or equivalent).
• Internal system error (which indicates a bug in Maple itself).
• ASSERT or return type or local variable type assertion failure (assertion failures are not trappable because they indicate a coding error, not an algorithmic failure).
Stack overflow (when that happens, generally stack space is insufficient for doing anything like running cleanup code). This includes the "too many levels of recursion" exception at the Maple language level.

I think...

Answered: tomleslie 13876

November 18 2022

1 0

you are going to have to upload an actual worksheet (using the big green up-arrow in the Mapleprimes toolbar) illustrating the problem, because the attached just works! (although for some reason it won't display inline on this site)

Download oddInt.mw

Well...

Answered: tomleslie 13876

November 18 2022

2 0

simply by using 'assuming x>0' on the evaluation of 'f', the attached (which for some reason won't display inline here) computes values for N=8 in about 20 seconds

badInt.mw

I assume...

Answered: tomleslie 13876

November 18 2022

0 0

you want something like the animation shown in the attached?

>

  restart:
  with(plots):
  with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':
  a := 7:
  b := 3:
  ellipse(p, x^2/a^2 + y^2/b^2 - 1, [x, y]):
  Fig := proc(t)
              local M, l1,L1,L2, q, m, c , tex;
              global a, b, p;
              point(M, a*cos(t), b*sin(t));
              if   VerticalCoord(M)=0
              then line(l1, x=HorizontalCoord(M));
              else m:=eval[recurse]( rhs
                                     ( isolate
                                       ( diff
                                         ( eval
                                           ( Equation(p),
                                             y=y(x)
                                           ),
                                           x
                                         ),
                                         diff(y(x),x)
                                       )
                                     ),
                                     [ y(x)=y,
                                       x=HorizontalCoord(M),
                                       y=VerticalCoord(M)
                                     ]
                                   );
                   c:=VerticalCoord(M)-m*HorizontalCoord(M);
                   line(l1, y=m*x+c);
              fi;
              reflection(q, p, l1);
              line(L1, [point(P1, evalf~(coordinates(foci(q)[1]))),
                            point(P2, evalf~(coordinates(foci(q)[2])))]):
              line(L2, [point(P3, evalf~(coordinates(foci(p)[1]))),
                            point(P4, evalf~(coordinates(foci(p)[2])))]):
              segment(P1P3,[P1,P3]):
              circle(cir1,[P3,2*a]):
              segment(P2P4,[P2,P4]):
              circle(cir1,[P3,2*a]):
              circle(cir2,[P4,2*a]):
              display
              ( [ textplot( [ [ HorizontalCoord(P2), VerticalCoord(P2), "F1"],
                               [ HorizontalCoord(P4), VerticalCoord(P4), "F2"]
                             ],
                             font = [times, 10], align = {below, left}
                           ),
                   draw( [ p(color = blue),P1P3(color=magenta,linestyle=3),
                           q(color = blue),P2P4(color=magenta,linestyle=3),
                           l1(color = black),L1(color = black),cir1(color=magenta,linestyle=3),
                           cir2(color=magenta,linestyle=3),
                           M(color = blue, symbol = solidcircle, symbolsize = 16),
                           P1(color = red, symbol = solidcircle, symbolsize = 16),
                           P2(color = red, symbol = solidcircle, symbolsize = 16),
                           P3(color = red, symbol = solidcircle, symbolsize = 16),
                           P4(color = red, symbol = solidcircle, symbolsize = 16)
                         ],
                         axes = none,
                         scaling = constrained
                       )
               ]
             );
          end proc:
#
# Change nFig to type float
#
   nFig := 180.0:
#
# Make sure loop index is a float, and change end point
# on loop to nFig-1.0
#
   Figs := seq(Fig(2*Pi*i/nFig), i = 0.0 .. nFig-1.0):
   display(Figs, insequence = true);

>

Download ellAgain.mw

Just for fun...

Answered: tomleslie 13876

November 17 2022

1 1

This is one of "the other ways" Acer hinted at

>

restart;

>

  with(SignalProcessing):
  with(plots):
  omega:= 10*2*Pi:
  X:= 0.01:
  Y:= 0.5:
  F:= z-> local n:
          add(n/5*cos(n*omega*z),n=1..5):
  M:= Matrix( round(Y/X)+1,
              2,
              (i,j)-> `if`( j=1,
                            (i-1)*X,
                            F((i-1)*X)
                          )
            ):
  pointplot(M, connect=true);
  DFT(M[..,2]);

(1)

>

Download DFT.mw

A random guess...

Answered: tomleslie 13876

November 17 2022

0 0

Since you do not specify where you want to

"rite lettres F1, F2"

For the attached, I just guessed that you might want to label the foci of the rotating ellipse.

If this is not what you want then you have to specify exactly what you want written and where!

>

  restart:
  with(geometry):
  _EnvHorizontalName := 'x':
  _EnvVerticalName := 'y':
  a := 7:
  b := 3:
  ellipse(p, x^2/a^2 + y^2/b^2 - 1, [x, y]):
  Fig := proc(t)
              local M, l1,L1,L2, q, m, c;
              global a, b, p;
              uses plots:
              point(M, a*cos(t), b*sin(t));
              if   VerticalCoord(M)=0
              then line(l1, x=HorizontalCoord(M));
              else m:=eval[recurse]( rhs
                                     ( isolate
                                       ( diff
                                         ( eval
                                           ( Equation(p),
                                             y=y(x)
                                           ),
                                           x
                                         ),
                                         diff(y(x),x)
                                       )
                                     ),
                                     [ y(x)=y,
                                       x=HorizontalCoord(M),
                                       y=VerticalCoord(M)
                                     ]
                                   );
                   c:=VerticalCoord(M)-m*HorizontalCoord(M);
                   line(l1, y=m*x+c);
              fi;
              reflection(q, p, l1);
              line(L1, [point(P1, evalf~(coordinates(foci(q)[1]))),
                            point(P2, evalf~(coordinates(foci(q)[2])))]):
              line(L2, [point(P3, evalf~(coordinates(foci(p)[1]))),
                            point(P4, evalf~(coordinates(foci(p)[2])))]):
              segment(P1P3,[P1,P3]): circle(cir1,[P3,2*a]):segment(P2P4,[P2,P4]):
              circle(cir1,[P3,2*a]): circle(cir2,[P4,2*a]):
              display( [ textplot
                         ( [ [ coordinates(foci(q)[1])[],
                               "F1",
                               align={above, right}
                             ],
                             [ coordinates(foci(q)[2])[],
                               "F2",
                               align={above, right}
                             ]
                           ]
                         ),
                         draw
                         ( [ p(color = blue),P1P3(color=magenta,linestyle=3),
                             q(color = blue),P2P4(color=magenta,linestyle=3),
                             l1(color = black),L1(color = black),cir1(color=magenta,linestyle=3),
                             cir2(color=magenta,linestyle=3),
                             M(color = blue, symbol = solidcircle, symbolsize = 16),
                             P1(color = red, symbol = solidcircle, symbolsize = 16),
                             P2(color = red, symbol = solidcircle, symbolsize = 16),
                             P3(color = red, symbol = solidcircle, symbolsize = 16),
                             P4(color = red, symbol = solidcircle, symbolsize = 16)
                           ],
                           axes = none,
                           scaling = constrained
                         )
                       ]
                     );
          end proc:
#
# Change nFig to type float
#
   nFig := 120.0:
#
# Make sure loop index is a float, and change end point
# on loop to nFig-1.0
#
   Figs := seq(Fig(2*Pi*i/nFig), i = 0.0 .. nFig-1.0):
   plots:-display(Figs, insequence = true);

>

Download ell.mw

Well...

Answered: tomleslie 13876

November 16 2022

0 2

as the help (rather unhelpfully) says (emphasis added)

Integrals appearing in answers returned by dsolve are expressed using the inert Int and Intat (not int or intat). These integrals appear when int is not able to calculate them or when it appears to be convenient not to evaluate them during the solving process. You can request the evaluation of these integrals using the value command.

So just add a value() command to your first example, as in the attached.

Download odeVal.mw

So many problems...

Answered: tomleslie 13876

November 14 2022

0 2

The first of the attached files is the code you supplied with added comments just to illustrate why there are o mant errors and ambiguities that it has no chance of ever working

Download OPevents2.mw

Your original question states

Basically there is a very simply system of two differential equations and I would simply like to see an event triggered when the dependent variable is above 2.

Well you don't state which independent variable, bu the following code will execute, and the dsolve() soltion proces will terminate whenever either dependent variable exceeds 2.

  restart;
  fcns := {y(x), z(x)}:
  p := dsolve( { diff(y(x), x) = 5*z(x),
                 diff(z(x), x) = 4*y(x),
                 y(0) = 0,
                 z(0) = 1
               },
               fcns,
               numeric,
               events = [ [y(x) - 2, halt],
                          [z(x) - 2, halt]
                        ]
             ):
  plots:-odeplot( p,
                  [ [x, y(x)],
                    [x, z(x)]
                  ],
                  x=0..0.3,
                  color=[red, blue]
                );

Download odeEvents2.mw

Well...

Answered: tomleslie 13876

November 13 2022

4 0

In order to plot the point F2 - you will have defined F2! You did define F1 twice, so I'm guessing this is a simple typo

You also need plots[pointplot]()

Both problems fixed in the attached, along with a "neater" way to do it, using the geometry() package

>

restart:
_local(D):
f := (x, y) -> 3*x^2 - 3*y*x + 6*y^2 - 6*x + 7*y - 9:
coeffs(f(x, y)):
A, B, C, D, E, F := %:
theta := 1/2*arctan(B/(A - C)):
solve({-2*A*xc - B*yc = D, -B*xc - 2*C*yc = E}):
assign(%):
x := xcan*cos(theta) - ycan*sin(theta) + xc:
y := xcan*sin(theta) + ycan*cos(theta) + yc:
Eq := simplify(expand(f(x, y))):
xcan^2/simplify(sqrt(-tcoeff(Eq)/coeff(Eq, xcan^2)))^`2` + ycan^2/simplify(sqrt(-tcoeff(Eq)/coeff(Eq, ycan^2)))^`2` = 1:
a := sqrt(-tcoeff(Eq)/coeff(Eq, xcan^2)):
b := sqrt(-tcoeff(Eq)/coeff(Eq, ycan^2)):
c := sqrt(a^2 - b^2):
F1 := [xc + c*cos(theta), yc + c*sin(theta)]:
evalf(%):
F2 := [xc - c*cos(theta), yc - c*sin(theta)]:
evalf(%):
Points := plots[pointplot]([F1[], F2[]], symbol = solidcircle, color = [red], symbolsize = 6):
xcan := plot(yc + tan(theta)*('x' - xc), 'x' = -2 .. 3.5, color = black):
ycan := plot(yc - ('x' - xc)/tan(theta), 'x' = 0.1 .. 1.5, color = black):
Ellipse := plots[implicitplot](f('x', 'y'), 'x' = -2 .. 3.5, 'y' = -2 .. 1.5, color = red, thickness = 2, gridrefine = 5):
labels := plots[textplot]([[0.4, 1.3, "ycan"], [3.2, 0.75, "xcan"]], font = [TIMES, ROMAN, 14]):
plots[display](xcan, ycan, Points, Ellipse, labels, scaling = constrained);

>

  restart:
  with(plots):
  with(geometry):
_EnvHorizontalName := 'x':
_EnvVerticalName := 'y':

  ellipse(el, 3*x^2 - 3*y*x + 6*y^2 - 6*x + 7*y - 9):
  evalf~(coordinates~(foci(el)));

  line( L1,
        [ point(P1, evalf~(coordinates(foci(el)[1]))),
          point(P2, evalf~(coordinates(foci(el)[2])))
        ]
      ):
  PerpendicularLine( PL1,
                     point( P3,
                            coordinates( center(el))
                          ),
                     L1
                   ):

  display( [ textplot( [ [rhs~(solve( [Equation(el), Equation(L1)], [x,y])[1])[], "major"],
                         [rhs~(solve( [Equation(el), Equation(PL1)], [x,y])[1])[], "minor"]
                       ],
                       'align'={'right'}
                     ),

             draw( [ el(color=red),
                     L1(color=black),
                     PL1(color=black),
                     P1(color=red, symbol=solidcircle, symbolsize=16),
                     P2(color=red, symbol=solidcircle, symbolsize=16)
                   ],
                   axes=normal
                 )
           ]
         );

>

Download ellipse.mw

Your fundamental problem...

Answered: tomleslie 13876

November 10 2022

0 0

is that when you *think* you are invoking the procedure, you actually have the code

TEST*(A, B, TOL, maxIter)

Notice the multiplication sign?!!!! This occurs because you have a space between TEST and (A, B, TOL, maxIter). So the 2D input parser regards this as "implicit multiplication", so inserts a multiplication symbol, ie '*'. Aah the pleasures of using 2D Math input!!

There are multiple instances of the same problem within the procedure itself, usually associated with the print() command, where you have typed print () rather than print() Again the extra space is the problem. It is tricky to find/fix all of these because (obviously) the standard find/replace capability does not work for 2D math input

As a side note, I notice a lot of return commands in your procedure which do not actually return anything. Is this deliberate?

I think I have fixed the problem with "extra" spaces in the attached, although I cannot be sure. The procedure now "executes" - but the problem is that it never seems to stop - I gave it ten minutes on my machine and it was still running! You could try the attached and see if it ever stops executing!!

Any further diagnosis probably requires firing up the debugger which I don't have time for now - past my bedtime!