Carl Love

Carl Love

28095 Reputation

25 Badges

13 years, 100 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are answers submitted by Carl Love

How about this?

degrees:= m-> map(x-> x=degree(m,x), indets(m)):
degrees(x*y^2);

Here's an object-oriented implementation of circular lists that allows for all of the multitude of types of indexing that ordinary lists allow (including unevaluated indices) and has as transparent an interface as I could manage. That is, I tried to make everything look like ordinary lists.


An object-oriented implementation of circular lists.

restart:


module CircularList()
option
     `Author: Carl Love, 2015-Jul-29`,
     object
;
export
     List, #the underlying list

     ModuleType::static:= proc(CL,t1,t2,$)
          nargs=3 implies CL:-List::list(t2)
     end proc,

     ModuleApply::static:= proc()
     local CL:= Object(CircularList);
          CL:-List:= [args];
          CL
     end proc,

     #Make circular lists prettyprint just like lists.
     ModulePrint::static:= CL-> CL:-List,

     `?[]`::static:= proc(CL,K)     
     local L:= CL:-List, N:= nops(L), k:= K[], n, _k;

           #Allow all types of index specifications that ordinary lists do.
           if not K::{
                          identical([]),
                          [{
                               list,
                               integer,
                               And(algebraic, Not(complexcons)),
                               range({integer, And(algebraic, Not(complexcons))})                             
                          }]
                     }

           #Return same error message as ordinary lists.
           then error "invalid subscript selector"

           #NULL index extracts the underlying sequence.
           elif K=[] then L[]
           
           #A list of indices returns a sub-CircularList.
           elif k::list then ModuleApply(map2(thisproc, CL, `[]`~(k))[])

           #A range of integers returns a sub-CircularList, but the range must be
           #a valid index for the underlying list.
           elif k::range(integer) then ModuleApply(L[k][])

           #If the CircularList is empty, return NULL.
           elif N=0 then

           #For plain integers, use modular arithmetic.
           elif k::integer then
                n:= modp(k,N);             
                if n=0 then n:= N end if;
                L[n]

           #The index has an unevaluated component.
           else             
                subs(_k= k, 'CL[_k]')
           end if
     end proc
;
end module:    
      

Examples

L1:= CircularList(x1,x2,x3,x4);

thismodule

L2:= CircularList(y1,y2,y3,y4);

thismodule

L1[7];

x3

L2[0];

y4

Just like ordinary lists, unevaluated indices are allowed.

x:= L1[n^2];

(thismodule)[object, `Author: Carl Love, 2015-Jul-29`]

eval(x, n= 3);

x1

type(L1, CircularList);

true

type(L1,'CircularList'(name));

true

type(L1,'CircularList'(posint));

false

Empty circular lists are allowed.

L3:= CircularList();

thismodule

They return NULL for any single specific index.

x:= L3[3];

"x:="

Garbage indices are rejected with the same message as for ordinary lists.

L1[.3];

Error, (in CircularList:-?[]) invalid subscript selector

 

A NULL index extracts the underlying sequence.

L1[];

x1, x2, x3, x4

Ranges of indices must be such that they work for the underlying list (no modular arithmetic is done). They return a sub-CircularList

L1[2..4];

thismodule

Lists of indices may be modular. They also return a sub-CircularList.

L2[[7,3,1]];

thismodule

Ranges with unevaluated parts are allowed.

L2[1..m];

(thismodule)[object, `Author: Carl Love, 2015-Jul-29`]

And they return a sub-CircularList when evaluated.

eval(%, m=3);

thismodule


Download CircularLists.mw

Replace solve with eliminate. (The arguments can stay the same.) Remove the assign command. The solution returned by eliminate comes in two parts, each a set. The first set is the solutions for variables that you specified in the second argument. Inspect it very carefully. Is the form what you expected? Make sure that you spelled the variable names the same in all use instances. Maple is case sensitive. I'm somewhat skeptical that you spelled them all correctly because the solution looks weird to me. The second set returned by eliminate is a set of expressions which must be equal to 0 in order for the first set to actually be a solution.

It is only slow if you try to display the results in the Standard GUI. There is no reason to display 17550 3x3 listlists. Just end your statement with a colon rather than a semicolon (to suppress display) and the allstructs will run in 0.1 seconds.

There are very small discontinuities in the piecewise function due to roundoff errors during its creation. The solution is to create the piecewise function at a higher Digits setting than is used during the fsolve. For example, setting Digits:= 20 before defining the function and then changing to Digits:= 15 before doing the fsolve works.

What have you tried?? The obvious thing to try is factor, which will tell you immediately that the polynomial can't be factored.

Several other points:

  • The polynomial isn't huge by Maple standards.
  • I don't think that you know what it means to know something a priori. It means to know something because it has been proven in the mathematical (or pure logic) sense of proof. That seems to be precisely what you don't have in this case. To know something by experience or by observation is to know it a posteriori.
  • There's no need to attempt a partial factorization in this case. The command factor will tell you in a blink of an eye that the polynomial can't be factored at all.
  • Please post your expressions with explicit multiplication signs--- * ---so that they can be copy-and-pasted directly into Maple.

In Maple 18, if you omit the boundary conditions (Venkat's suggestion) and cancel the extraneous exponential factor in the first ODE, then dsolve will return a solution. In the worksheet below, the single line of 1-D input is my addition to cancel the factor. 

Note that dsolve's solution is explicit although at first glance it may seem implicit.


restart

n := 0:

ode1 := diff(exp(-beta*theta(y))*(diff(u(y), y)), y) = 0;

-beta*(diff(theta(y), y))*exp(-beta*theta(y))*(diff(u(y), y))+exp(-beta*theta(y))*(diff(diff(u(y), y), y)) = 0

ode1:= simplify(ode1*exp(beta*theta(y)));

-(diff(u(y), y))*(diff(theta(y), y))*beta+diff(diff(u(y), y), y) = 0

ode2 := diff((1+theta(y))^n*(diff(theta(y), y)), y)+Br*exp(-beta*theta(y))*(diff(u(y), y))^2 = 0;

diff(diff(theta(y), y), y)+Br*exp(-beta*theta(y))*(diff(u(y), y))^2 = 0

sol := dsolve([ode1, ode2])

[{theta(y) = ln(1+tan((1/2)*(_C2*beta)^(1/2)*(_C3+y)*2^(1/2))^2)/beta+_C4}, {u(y) = Int((-Br*exp(-beta*theta(y))*(diff(diff(theta(y), y), y)))^(1/2)/(Br*exp(-beta*theta(y))), y)+_C1, u(y) = Int(-(-Br*exp(-beta*theta(y))*(diff(diff(theta(y), y), y)))^(1/2)/(Br*exp(-beta*theta(y))), y)+_C1}]


Download Parallel_flow.mw

 

It's all about the order of evaluation. One way around the problem is to exploit the "special evaluation rules" of seq by simply replacing eval with seq:

seq(g(x), n= 3);

You could just as well use add or mul. These will work in Maple versions older than eval[recurse].

member(4, A);

 

Suppose that your list of equations is copied from Matlab as a single string, which I'll call MS (Matlab String), like this

MS:= "C_p_e = C_state/C_c;
C_p_f = I_state/I_i;
R_p_e = R_r*C_p_f;
I_p_e = (Se_p_e-R_p_e)-C_p_e;
";

The presence of newline characters in the string is insigificant---Maple ignores them. The semicolons, however, are necessary. Then do

eq:= parse~(StringTools:-Split(MS, ";"));

Now eq is a list of all the equations, so eq[1] is C_p_e = C_state/C_c, etc. You just need to use subs once, on all of eq. It's more efficient than calling allsubs repeatedly.

eq_a:= subs(
     C_p_e=fs(t), C_p_f=v(t), R_p_e=fd(t), I_p_e=fm(t), C_c=1/k, I_i=m, 
     R_r=c, Se_p_e=fe(t), I_state=int(fm(t),t), C_state=int(v(t),t),
     eq
);

The solutions are not necessarily real in your example. If the solutions were necessarily real, you could use syntax almost identical to the Mathematica:

a:= INTERVAL(2..4):
b:= INTERVAL(10..11):
c:= INTERVAL(1..3):
evalr~([solve(a*x^2+b*x+c)]);

If you want to do complex range arithmetic with your original example, you could do this:

a:= INTERVAL(2,0,4,0):
b:= INTERVAL(1,0,3,0):
c:= INTERVAL(1,0,3,0):
evalrC~([solve(a*x^2+b*x+c)]);

However, the results are not perfect, because they are in floating point. Note how a complex interval is specified with four numbers. See ?evalrC and ?evalr.

simplify(x1=a-y1-d*y2, {a-y2-d*y1= x2, 1-d^2= b, a-a*d= c});

Note that the three substitution equations must be written with the expressions on the left and the single variables on the right.

Don't use map when applying expand or simplify to a single expression. Use map to distribute these operations over a Matrix, but not on an single entry extracted from a Matrix.

There are also "tutors" for this in the Student package, but I think that it's instructive to see the components of a plot command.

 

restart:

 

Your first question: Plotting a function and one of its tangent lines.

f:= x-> 2*cos(x)-x^2:  a:= 2:  r:= 3:

 

Here I use a columnar layout to emphasize the parallelism between the plot components and their attributes.

plot(
                 [f(x),  D(f)(a)*(x-a)+f(a), [[a,f(a)]]], x= a-r..a+r,
      style=     [line,  line,               point     ],
      color=     [red,   blue,               green     ],
      linestyle= [solid, dot,                solid     ],
      symbolsize= 16, symbol= solidcircle,
      gridlines= false
);

Your second question: Plotting a rational function, its derivative, vertical asymptotes, and critical points.

f:= x-> (x^3-10*x^2-2*x+1)/(4*x^3+5*x+1);

proc (x) options operator, arrow; (x^3-10*x^2-2*x+1)/(4*x^3+5*x+1) end proc

Find the vertical asymptotes as the zeros of the denominator.

VA:= fsolve(denom(f(x))=0);

-.194145720502372

Find the critical points as the zeros of the numerator of the derivative.

CP:= fsolve(numer(D(f)(x))=0);

-1.47480661500595, 1.14215480870403

Get the horizontal domain, including all x-values of interest)...

(a,b):= (min,max)(VA,CP);

-1.47480661500595, 1.14215480870403

...and stretch it a little.

HR:= a-.3*(b-a) .. b+.3*(b-a):

 

And the same for the vertical range.

(a,b):= (min,max)(f~([CP]))[];  YR:= a-.3*(b-a) .. b+.3*(b-a);

-1.01333055349349, 1.09390932026410

-1.64550251562077 .. 1.72608128239138

plots:-display([
     #Use option discont for plots with vertical asymptotes:
     plot([f(x), D(f)(x)], x= HR, discont, legend= ["f", "f '"]),
     #Use parametric plots for vertical lines:
     map(va-> plot([va, t, t= YR], color= red, linestyle= dash), [VA])[],
     map(cp-> plot([[cp,f(cp)]], style= point, color= green), [CP])[]
     ],
     view= [DEFAULT, YR], gridlines= false,
     symbol= solidcircle, symbolsize= 16
);

 

 

Download RatFuncPlot.mw

 

In the future, you should show what you tried. Here's how to plot it.

f:= x-> 2*cos(x)-x^2:  a:= 2:  r:= 3:
plot([f(x), D(f)(a)*(x-a)+f(a)], x= a-r..a+r, color= [red,blue], linestyle= [solid,dot]);

I hope that you can generalize it from there.

First 254 255 256 257 258 259 260 Last Page 256 of 395