Carl Love

Carl Love

26827 Reputation

25 Badges

11 years, 281 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

Your matrix is severely rank deficient, a 57x57 with rank 44. The matrix is singular. I don't think it's a problem in your program; rather, it's a problem with your initial data. You could seek a least-squares solution to the linear system. However, I know nothing about this type of engineering problem, so I don't know if that's appropriate.

You should not use the package linalg. Use LinearAlgebra instead.

By the way, this is the most primitive Maple code of its size that I have ever seen---although it is neatly presented. It's very easy to make a mistake if you initialize huge matrices by indexing each element in its own assignment statement.

You wrote:

can this function be animated in 3d?

Sure. It's very easy:

S:= pdsolve({pde, bc}):
plots:-animate(plot3d, [eval(U, S), x= 0..1, t= 0..T], T= 0..10);

This particular function makes a rather boring animation though.

Why are you using Vectors instead of lists? You can use Vectors, but for every operation of this type, you'll need to convert to a list and then back to a Vector. And possibly you'll need to convert to and from a set also.

Once you have it in list form, there are several commands in ?ListTools to help with the things that you're asking about, specifically, ?ListTools,Categorize ?ListTools,Classify ?ListTools,FindRepetitions ?ListTools,MakeUnique and ?ListTools,SearchAll

You do it by assigning the output of fsolve to a variable (let's say Sol), and then using Sol as the second argument to eval with the first argument being the expression which contains the variables that fsolve solved for. (This same technique can be applied to any of the "solve" commands: dsolve, solve, rsolve, pdsolve, isolve, msolve, etc.)

Example:

Sys:= {'randpoly([x||(1..3)])+randpoly(x, degree= 0)' $ 3}; 
 /           3        3   2       3           3              3
{ 31 x1 x2 x3  - 27 x2  x3  - 2 x1  x2 - 31 x2  x3 + 25 x2 x3
 \                                                            

                            4        5       4                 2
   - 97 x2 x3 + 65, 68 x1 x2  - 29 x2  + 5 x2  x3 + 73 x1 x2 x3

             3           2            5        4   
   + 95 x1 x3  + 31 x1 x3  - 26, 91 x1  - 22 x1  x2

                3        2   2        4                   \
   + 51 x1 x2 x3  + 52 x1  x3  + 36 x3  + 18 x1 x2 x3 - 27 }
                                                          /
Sol:= fsolve(Sys, {x||(1..3)});
      {x1 = -0.194514327765446, x2 = -0.806303974671123,
    x3 = -0.998794548580666}

plot(eval(x1*t^2 + x2*t + x3, Sol), t= -2..2);

To add static frames to an animation, you first create the animation (using, for example, plots:-animate or plots:-display(..., insequence= true)), and then you use display again to merge the animation with a single static plot, which becomes a fixed part of every frame. This merged animation can be played as is, or it can be merged with another animation, once again with didplay. There's no limit to how many levels of merging with display that you can do. So, you can create animations that have parts that are static during only part of the animation.

So, here is your code with my modifications. It uses four levels of display, whereas your original had two levels:

restart: 
with(plots):
auto:= (x, y, C)-> plots:-pointplot([[y, x]], color= C):

p1:= animate(auto, [t, .17+0.55e-1*t, blue], t= .2 .. -2, frames= 10):
p2:= animate(auto, [t, .16+0.45e-1*t, red], t= .7 .. -2, frames= 20):
p3:= animate(auto, [t, .15+0.45e-1*t, green], t= 1.2 .. -2, frames= 30):
p4:= animate(auto, [t, .15+0.45e-1*t, black], t= 1.7 .. -2, frames= 40):

p5:= animate(auto, [t, .12+0.45e-1*t, blue], t= 2.2 .. .2, frames= 50):
lastframe5:= auto(eval([t, .12+0.45e-1*t, blue], t= .2)[]):

stoppt:= frame-> 2.7+(.7-2.7)*frame/60: #Split t-range for auto 6.
p61:= animate(auto, [t, .11+0.45e-1*t, red], t= 2.7 .. stoppt(50), frames= 50):
#Frames 51-60 for auto 6:
p62:= animate(auto, [t, .11+0.45e-1*t, red], t= stoppt(51) .. .7, frames= 10):
#Merge animations for auto 6:
p6:= display([p61, display([p62, lastframe5])], insequence):
lastframe6:= auto(eval([t, .11+0.45e-1*t, red], t= .7)[]):

firstframe7:= auto(eval([-.17+0.35e-1*t, t, green], t= -.2)[]):
p7:= animate(auto, [-.17+0.35e-1*t, t, green], t= -.2 .. 7, frames= 20):

firstframe8:= auto(eval([-.15+0.25e-1*t, t, black], t= -.7)[]):
p8:= animate(auto, [-.15+0.25e-1*t, t, black], t= -.7 .. 6, frames= 20):

A:= display([firstframe||(7..8), p||(1..6)]):
B:= display([p||(7..8), lastframe||(5..6)]):
display(
   [A, B], insequence,
   scaling= constrained, symbol= solidbox, symbolsize= 20,
   title= ``
);

The improved animation:

Your original animation:

I have a feeling that this can be handled by the events option to dsolve (see ?dsolve,events ), but I don't know the details. There are far too few examples on the help page and an overwhelming number of options described above the examples.

But here is a quick-and-dirty (linear search) procedure that will give you the last values plotted by odeplot (the last X value and the last Y value) before the first undefined values. (I'd be pleased if someone else could comment on a better way to do this.)

LastValues:= proc(P::specfunc(anything,PLOT))
     local M,X,Y,m,n;
     M:= op([1,1], P);
     X:= convert(M[.., 1], list);
     m:= ListTools:-Search(HFloat(undefined), X);
     Y:= convert(M[.., 2], list);
     n:= ListTools:-Search(HFloat(undefined), Y);
     `if`(m < 2, undefined, X[m-1]), `if`(n < 2, undefined, Y[n-1])
end proc;

Call the procedure as, for example, LastValues(plott[1]), and it should return two numbers, the first being the last X value and the second being the last Y value.

Please let me know if this does what you need.

dotprod(j,k) is 0 after you radnormal or simplify it. There is nothing wrong with your code.

You should not use linalg anymore. Use LinearAlgbera. It is available in Maple 8.

I can't tell if you want to create a surface, or if you simply want the sequence of curves in a 3D plot. If you just want the curves, it requires merely a trivial modification to your code:

  1. Change the curve that you plot in the odeplot command from [X(t), Y(t)] to [ j, Y(t), X(t)].
  2. Change the display command to display(convert(plott, list)).
  3. (Optional) Include an option such as thickness= 5 in the odeplot to improve visibility.

If you do want to create a surface, let me know.

 

[I converted your Post into a Question.-- Carl Love as a moderator]

Please upload a worksheet and/or post your equations in plaintext so that we do not need to retype them into Maple. Thanks.  The solution to your problem is probably to simply include a range in the fsolve command. See ?fsolve,details

fsolve({f, g}, {x= 0..Pi/2, y= 0..Pi/2});

Here's how to do it with a continuous transformation to your existing color function, which is presumed to return a value between 0 and 1 (the HUE color scale). Keeping it continuous is very very nice when you want colors to represent  numeric values. Let's say your existing color function is C, and your coordinate functions for a parametrized surface are Fx, Fy, Fz.

Gamma:= 1.15:
plot3d(
     [Fx, Fy, Fz],  a..b, c..d,
     color=  [
          (x,y)-> (1-C(x,y))^Gamma/3, #Hue
          (x,y)-> 1-C(x,y)/4,         #Saturation
          (x,y)-> 1-C(x,y)/7,         #Value
          colortype= HSV
     ],
     lightmodel= NONE,
     style= patchnogrid     
);

There are several parameters that can be adjusted; I've chosen some of them by my personal taste for color .

  • Gamma controls the evenness of the distribution between red and green. I gave this one a name because this is a well-known concept (see the Wikipedia article "Gamma correction").
  • The 3 in the Hue selects the fraction (1/3 in this case) of the full color spectrum that you want. If you want green to red, it will need to be pretty close to 3.
  • The Hue value is subtracted from 1 to make the scale go green to red rather than red to green.
  • The 4 in the Saturation controls (to some extent) how "light" the light-green is.
  • The 7 in the Value controls (to some extent) how dark the dark-red is (lower values will make it darker).
  • lightmodel= NONE is used so that the colors will not change due to shadows when the plot is rotated.

Here is the resulting sub-spectrum:

Brian wrote:

the above solution violates my assumption phi should be between -1.57 and +1.57. why?

Markiyan's Answer addresses well how to get the results. Regarding why you had the problem: solve generally ignores assumptions (see ?solve,details), sometimes without warning. Instead of using assumptions, you can include inequalities with the equations that you pass to solve, as Markiyan did.

You wrote:

so I dont actually need to asign as I did with xs??

Like virtually all commands in Maple, you do need to assign the results of dsolve in order to use them later.

You wrote:

in ordinary circumstances plot(eval(x(t),t=0..3) would be sufficient?

No, the command is plot(eval(x(t), xs), t= 0..3). It requires the xs. Note that you cannot directly plot the results of a dsolve(..., numeric); for the numeric case, use plots:-odeplot instead.

You wrote:

could I similarly use eval(xs,x(t)) ??

No, the command is eval(x(t), xs). The second argument to this type of eval must be an equation, or a list or set of equations--- which is what xs is in the case of a non-numeric dsolve.

You wrote:

ode3 := diff(y(t), t, t) = y(t)*sqrt((1-y(t))^2-1)/(1-y(t))

ics3 := y(0) = -2, (D(y))(0) = 0

ys := dsolve({ics3, ode3}, numeric)

I can ode plot this as the first one.

But my next step was to parametrically plot both solutions like

plots:-odeplot([xs(t), ys(t), t = 0 .. 3])

plots:-odeplot([xs, ys, t = 0 .. 3])

plots:-odeplot([x(t), y(t), t = 0 .. 3])

Unfortunately, each invocation of plots:-odeplot can only be used with the results of a single dsolve(..., numeric). But there are two ways around this problem. The first is to apply dsolve to the ODEs together as a system:

Sol:= dsolve({ode3, ics3, ode, ics}, numeric):
plots:-odeplot(Sol, [[t,x(t)], [t,y(t)]], t= 0..3);

The second is to combine two plots with plots:-display (which can combine any number of plots):

plots:-display([
     plots:-odeplot(xs, t= 0..3),
     plots:-odeplot(ys, t= 0..3)
]);

The second method gives you more flexibility.

You wrote:

how can one plot:

plots:-odeplot(f(xs), t= 0..3);

with some function f

Like this:

plots:-odeplot(xs, [t, f(x(t))], t= 0..3);

The command that you need is

B:= copy(A);

See ?copy for some details why. A few types of "container" structures in Maple are considered "mutable": tables, rtables (which includes Arrays, Matrices, and Vectors), procedures, and modules (which includes Records). Mutable objects need to be fully copied in order for an assignment to one to not affect the other. Note that most structures in Maple are not mutable even if they are containers: lists, sets, sequences.

 

As you have written it, it's an ODE, not a PDE, because all the derivatives are taken with respect to (wrt) x. Perhaps you meant for the derivative on the right side of the equals sign to be wrt t instead? Since it's a second-order ODE wrt x, two boundary conditions are enough, hence the error message.

The result returned by your dsolvexs, is of the form x(t) = ....  The plot command is complaining because you passed it an equation. The actual solution is the right-hand side of the equation. There are many ways  to isolate the right-habd side. One way that's easy to remember is rhs(xs). The way I generally prefer is eval(x(t), xs). Under ordinary circumstances you could give the plot command

plot(eval(x(t), xs), t= 0..3);

But these are not ordinary circumstances. The solution, although it appears short, is so complicated that Maple is having a lot of trouble evaluating it numerically. What we need to do is get a numeric solution from dsolve and then use plots:-odeplot to plot it:

xs:= dsolve({ics, ode}, numeric);
plots:-odeplot(xs, t= 0..3);

By the way, the right side of your ode has the useless coefficient 1^2. Did you intend for that to be something else?

First 356 357 358 359 360 361 362 Last Page 358 of 385