Carl Love

Carl Love

28085 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

It might be possible if the expressions are polynomials. See ?solve,parametric. The algorithm used by this command is extremely complicated, and results often take a long time to compute.

Using RealDomain seems to interfere with the ability of simplify to process assumptions. So, don't use it. I don't know exactly why it interferes, but I do know that RealDomain alters the definitions of both sqrt and simplify.

AFAIK, testeq is a probabalistic test that uses modular arithmetic. Because of the modularity, I don't think that it can ever process assumptions. You should replace testeq with is, which is a deterministic test that uses ordinary complex-number arithmetic and is the main command for deriving true/false conclusions from assumptions.

 

In the example below, is the summand, which I've represented as an m-ary function:

n_0:= 4:  m:= 3:
add(S(seq(i +~ 1)), i= Iterator:-Combination(n_0, m));

      S(1, 2, 3) + S(1, 2, 4) + S(1, 3, 4) + S(2, 3, 4)

The inner seq will always sequence the combinations in order.

What you referred to as "all monotonic lidyd [lists] of size m of [1..n_0]" is usually called all m-combinations of {1..n_0}.

An older-style command with the same result is

add(S(i[]), i= combinat:-choose(n_0, m));

The function osum could be written like

osum:= (S, k::(name= [nonnegint, nonnegint]))->
    local j;
    add(eval(S, lhs(k)= j+~1), j= Iterator:-Combination(rhs(k)[]))
:

And called like

osum(i[1]^2+i[2]^2+i[3]^2, i= [4,3]);
 

You could do this:

j:= piecewise(x < 0, x^2 + 1, 0 <= x, x - 1);
plot(j, x= -3..3, discont);

If I do as above, I can right click on the output of the first line, then "Open Context Panel", then "Plot Builder".

The error message is fairly clear: In order to get a numeric solution for x, you need to supply numeric values for hkpT[A], omega[A]. This is also mathematically clear. It seems unlikely that you're looking for a symbolic solution for x, because you've specified such a narrow numeric range for it.

I don't know what version of Maple you're using, but on mine it is already on the plot context menu. Right click on a plot for the context menu, then Symbol (7th item from the top in my version), then Symbolsize (bottom of the submenu), then enter the number.

The following procedure returns a root-bounding interval on which the original function is continuous, and it doesn't depend on the function or its numerator being polynomial.

p:= (x^3 - 13*x^2 + 55*x - 73)/(x - 1)
:
RootBoundwithDiscontinuity:= proc(f::algebraic, ab::name=range(realcons))
local A,k;
    for A in indets(plot(f, ab, discont), Matrix) do
        for k to upperbound(A,1)-1 do
            if A[k,2]*A[k+1,2] <= 0 then return [A[k,1],A[k+1,1]] fi
        od
    od;
    return
end proc
:
RootBoundwithDiscontinuity(p, x= -9..9); 
     [2.60815727428114, 2.64908352866142]

 

In addition to what Acer showed you, it's possible to combine any number of points, lines, and curves in a single plot command. Here's one with numerous bells & whistles.

L:= y = -3/2*x + 2; Pt:= [2,1]:
m:= -1/implicitdiff(L, y, x);
Lperp:= y = m*(x-Pt[1]) + Pt[2];
PtInt:= eval([x,y], solve({L, Lperp}));
plot(
    [rhs(L), rhs(Lperp), [PtInt], [Pt]], x= -5..5,
    style= [line$2, point$2], 
    symbolsize= 18, symbol= [solidcircle, soliddiamond],
    scaling= constrained, color= [red, "DarkGreen", black, red],
    view= [-5..5, -5..5], thickness= 0, size= [900$2],
    legend= [
        typeset(L, `\t`), 
        typeset(Lperp, ` `, `#mo(&perp;)`,`\t`),
        typeset(PtInt, `\t`),
        Pt
    ],
    labels= [x,y], 
    caption= "\nThe given information is in red.", 
    captionfont= [Times, 18],
    title= "Perpendicular Lines\n", titlefont= [Helvetica, Bold]
);

Here's a different way that uses neither complex arithmetic (like VV's) nor an idiosyncratic fractal-specific syntax (like Tom's). I have nothing against those methods; I just think that this real-arithmetic-only approach provides another valuable way of understanding what's happening.

N:= 6: #number of iterations
#Rotate point P by angle theta around point Q:
Rot:= theta-> (P,Q)-> 
    local R:= P-Q, s:= sin(theta), c:= cos(theta): 
    <R[1]*c - R[2]*s | R[1]*s + R[2]*c> + Q
:
Rot120:= Rot(2*Pi/3): #the only rotation that we need
#Select a point on line segment PQ:
BetweenPt:= (P,Q,r)-> P*(1-r) + Q*r:
#Replace line segment __ with _/\, i.e., put a notch in the middle:
Notch:= (P,Q)-> 
    local P1:= BetweenPt(P,Q,1/3); 
    (P, P1, Rot120(P,P1), BetweenPt(P,Q,2/3))
: 
Pts[1]:= <seq(<cos(2*Pi/3*k) | sin(2*Pi/3*k)>, k= 0..3)>:
for k to N do
    Pts[k+1]:= 
        <seq(Notch(Pts[k][j], Pts[k][j+1]), j= 1..3*4^(k-1)), <1|0>>
od:
plots:-display(
   seq(plot(Pts[k])$4, k= 1..N+1), insequence, 
   thickness= 0, scaling= constrained, axes= none
);

I made 4 changes to your worksheet. The change to the plot command is syntactically required.  The other 3 changes are just to increase the speed of numeric integration:

  1. Digits:= 15: #Increase from 10
  2. E:= (x,t)-> Int(  #Use uppercase I.
        exp((-esp*w^4+Disp*w^2+k)*t)*cos(w*(x+v*t))/Pi,
        w= 0..infinity,
        epsilon= 1e-7  #Added for speed.
    );
  3. #Same changes with larger epsilon:
    u:= (x,t)-> Int(E(x-xi, t)*f(xi), xi= 0..2e4, epsilon= 1e-5);
  4. plot(u(1500, t), t= 0..6e4, numpoints= 100);

This is the plot I got for numpoints= 10:

You'll need to wait several minutes for numpoints= 100.

In the worksheet below, the procedure represents your 2nd-order ODE as a 1st-order system in functional form. The procedure Heun specifies the iteration for what I believe your instructor means by "modified Euler". However, if this formula seems unfamiliar to you, then specify the formula that you're supposed to use, and this worksheet can be adapted by changing that 1 line.
 

restart
:

F:= (t,Y)-> <Y[2], -_c*Y[2] -_k*Y[1]>:

makeF:= (c,k)-> subs([_c= c, _k= k], eval(F))
:

Euler:= (F,t,Y,h)-> Y + h*F(t,Y):

Heun:= (F,t,Y,Yp,h)-> Y + h/2*(F(t,Y)+F(t+h, Yp))
:

Pred_Corr:= proc(F, Y0, t0, h, n, Pred, Corr, {digits::posint:= Digits})
local
    j, k, t:= t0, m:= numelems(Y0), Y:= copy(Y0), Yp:= copy(Y0),
    R:= Array(1..n+1, 1..2*m+1, [[t, seq(Y[j]$2, j= 1..m)]], datatype= sfloat)   
;
    UseHardwareFloats:= false;
    for k from 2 to n+1 do
        Yp:= evalf[digits](Pred(F,t,Y,h));
        Y:= evalf[digits](Corr(F,t,Y,Yp,h));
        t:= t+h;
        R[k,..]:= Vector[row]([t, seq([Yp[j], Y[j]][], j= 1..m)])
    od;
    R
end proc
:  

Case:= Record(t0= 0, y0= 1, yp0= 0, h= 0.02, n= 5, digits= 5):
Cases:= Record[Case]~(c=~ [0,2,2,4], k=~ [4,0,1,6.25]):
R:= table():

for C in Cases do
    R[c=C:-c, k=C:-k, h=C:-h]:= Pred_Corr(
        makeF(C:-c, C:-k), <C:-y0,C:-yp0>, C:-t0, C:-h, C:-n,
        Euler, Heun, digits= C:-digits
    )
od:

(P-> lhs(P)=
    <<t | `y predicted` | `y corrected` | `y' predicted` | `y' corrected`>, rhs(P)>
)~(<indices(R, pairs, indexorder)>);

Vector(4, {(1) = (c = 0, k = 4, h = 0.2e-1) = (Matrix(7, 5, {(1, 1) = t, (1, 2) = `y predicted`, (1, 3) = `y corrected`, (1, 4) = `y' predicted`, (1, 5) = `y' corrected`, (2, 1) = 0., (2, 2) = 1., (2, 3) = 1., (2, 4) = 0., (2, 5) = 0., (3, 1) = 0.2e-1, (3, 2) = 1., (3, 3) = .99920, (3, 4) = -0.8e-1, (3, 5) = -0.80000e-1, (4, 1) = 0.4e-1, (4, 2) = .99760, (4, 3) = .99680, (4, 4) = -.15994, (4, 5) = -.15987, (5, 1) = 0.6e-1, (5, 2) = .99360, (5, 3) = .99281, (5, 4) = -.23961, (5, 5) = -.23948, (6, 1) = 0.8e-1, (6, 2) = .98802, (6, 3) = .98723, (6, 4) = -.31890, (6, 5) = -.31872, (7, 1) = .10, (7, 2) = .98086, (7, 3) = .98007, (7, 4) = -.39770, (7, 5) = -.39744})), (2) = (c = 2, k = 0, h = 0.2e-1) = (Matrix(7, 5, {(1, 1) = t, (1, 2) = `y predicted`, (1, 3) = `y corrected`, (1, 4) = `y' predicted`, (1, 5) = `y' corrected`, (2, 1) = 0., (2, 2) = 1., (2, 3) = 1., (2, 4) = 0., (2, 5) = 0., (3, 1) = 0.2e-1, (3, 2) = 1., (3, 3) = 1., (3, 4) = 0., (3, 5) = 0., (4, 1) = 0.4e-1, (4, 2) = 1., (4, 3) = 1., (4, 4) = 0., (4, 5) = 0., (5, 1) = 0.6e-1, (5, 2) = 1., (5, 3) = 1., (5, 4) = 0., (5, 5) = 0., (6, 1) = 0.8e-1, (6, 2) = 1., (6, 3) = 1., (6, 4) = 0., (6, 5) = 0., (7, 1) = .10, (7, 2) = 1., (7, 3) = 1., (7, 4) = 0., (7, 5) = 0.})), (3) = (c = 2, k = 1, h = 0.2e-1) = (Matrix(7, 5, {(1, 1) = t, (1, 2) = `y predicted`, (1, 3) = `y corrected`, (1, 4) = `y' predicted`, (1, 5) = `y' corrected`, (2, 1) = 0., (2, 2) = 1., (2, 3) = 1., (2, 4) = 0., (2, 5) = 0., (3, 1) = 0.2e-1, (3, 2) = 1., (3, 3) = .99980, (3, 4) = -0.2e-1, (3, 5) = -0.19600e-1, (4, 1) = 0.4e-1, (4, 2) = .99941, (4, 3) = .99922, (4, 4) = -0.38812e-1, (4, 5) = -0.38424e-1, (5, 1) = 0.6e-1, (5, 2) = .99845, (5, 3) = .99827, (5, 4) = -0.56871e-1, (5, 5) = -0.56495e-1, (6, 1) = 0.8e-1, (6, 2) = .99714, (6, 3) = .99696, (6, 4) = -0.74201e-1, (6, 5) = -0.73835e-1, (7, 1) = .10, (7, 2) = .99548, (7, 3) = .99531, (7, 4) = -0.90821e-1, (7, 5) = -0.90466e-1})), (4) = (c = 4, k = 6.25, h = 0.2e-1) = (Matrix(7, 5, {(1, 1) = t, (1, 2) = `y predicted`, (1, 3) = `y corrected`, (1, 4) = `y' predicted`, (1, 5) = `y' corrected`, (2, 1) = 0., (2, 2) = 1., (2, 3) = 1., (2, 4) = 0., (2, 5) = 0., (3, 1) = 0.2e-1, (3, 2) = 1., (3, 3) = .99875, (3, 4) = -.1250, (3, 5) = -.12000, (4, 1) = 0.4e-1, (4, 2) = .99635, (4, 3) = .99520, (4, 4) = -.23524, (4, 5) = -.23048, (5, 1) = 0.6e-1, (5, 2) = .99059, (5, 3) = .98953, (5, 4) = -.33644, (5, 5) = -.33192, (6, 1) = 0.8e-1, (6, 2) = .98289, (6, 3) = .98192, (6, 4) = -.42906, (6, 5) = -.42476, (7, 1) = .10, (7, 2) = .97342, (7, 3) = .97254, (7, 4) = -.51352, (7, 5) = -.50944}))})

#Compare "corrected" columns with stock Maple command:
for C in Cases do
    (c=C:-c, k=C:-k, h=C:-h) = dsolve(
        {
            diff(y(t),t$2) + C:-c*diff(y(t),t) + C:-k*y(t),
            y(0) = C:-y0, D(y)(0) = C:-yp0
        },
        numeric, method= classical[heunform], stepsize= C:-h,
        output= Array(C:-t0 +~ C:-h*~[$0..C:-n])
    )
od;

(c = 0, k = 4, h = 0.2e-1) = (Vector(2, {(1) = Vector[row](3, {(1) = t, (2) = y(t), (3) = diff(y(t), t)}), (2) = Matrix(6, 3, {(1, 1) = 0., (1, 2) = 1., (1, 3) = 0., (2, 1) = 0.2e-1, (2, 2) = .9992, (2, 3) = -0.8e-1, (3, 1) = 0.4e-1, (3, 2) = .99680064, (3, 3) = -.15987200000000001, (4, 1) = 0.6e-1, (4, 2) = .992805759488, (4, 3) = -.2394881536, (5, 1) = 0.8e-1, (5, 2) = .9872217518084095, (5, 3) = -.3187210238361601, (6, 1) = .10, (6, 2) = .9800575539302396, (6, 3) = -.3974437871617639})}))

(c = 2, k = 0, h = 0.2e-1) = (Vector(2, {(1) = Vector[row](3, {(1) = t, (2) = y(t), (3) = diff(y(t), t)}), (2) = Matrix(6, 3, {(1, 1) = 0., (1, 2) = 1., (1, 3) = 0., (2, 1) = 0.2e-1, (2, 2) = 1.0, (2, 3) = 0., (3, 1) = 0.4e-1, (3, 2) = 1.0, (3, 3) = 0., (4, 1) = 0.6e-1, (4, 2) = 1.0, (4, 3) = 0., (5, 1) = 0.8e-1, (5, 2) = 1.0, (5, 3) = 0., (6, 1) = .10, (6, 2) = 1.0, (6, 3) = 0.})}))

(c = 2, k = 1, h = 0.2e-1) = (Vector(2, {(1) = Vector[row](3, {(1) = t, (2) = y(t), (3) = diff(y(t), t)}), (2) = Matrix(6, 3, {(1, 1) = 0., (1, 2) = 1., (1, 3) = 0., (2, 1) = 0.2e-1, (2, 2) = .9998, (2, 3) = -0.196e-1, (3, 1) = 0.4e-1, (3, 2) = .9992158800000001, (3, 3) = -0.3842384e-1, (4, 1) = 0.6e-1, (4, 2) = .9982629295600001, (4, 3) = -0.56494571952e-1, (5, 1) = 0.8e-1, (5, 2) = .9969559833638288, (5, 3) = -0.738346392364672e-1, (6, 1) = .10, (6, 2) = .9953094332381214, (6, 3) = -0.9046589172448144e-1})}))

(c = 4, k = 6.25, h = 0.2e-1) = Matrix(%id = 18446746390064941646)

 


 

Download PredictorCorrector.mw

If the integration is in rectangular coordinates, as it is in the example that you present, then a plot of the region of integration can be automatically generated from the inert form of the integral, like this:

J:= VectorCalculus:-int(x*y, [x,y]= Triangle(<0,0>, <1,0>, <0,1>), inert):
plot3d(0, op([1,2], J), op(2,J), orientation= [180,0,180]);

Then the integral can be computed by

value(J);

There's no need to create a 60,000 point spline! Just plot the points:

X:= ExcelTools:-Import("E:/ct.xls", "Sheet1"):
Y:= ExcelTools:-Import("E:/ct.xls", "Sheet2"):
plot(
    <Y | X[..,2]>, 
    labels= ['t', 'concentration'], labeldirections=[horizontal, vertical]
);

Both parts a) and c) can be answered with the same vector cross product. This can be accessed by loading the LinearAlgebra package and using the &x operator. It's as simple as this:

restart:
with(LinearAlgebra):
V:= <x, y, f(x,y)>: P:= <5,10,0>: Q:= <-1,7,1>: R:= <2,1,4>:
#Order doesn't matter as long as all points are used:
solve((nv:= (P-Q) &x (R-Q)).(P-V), {f(x,y)})[]; #linear function
Norm(nv,2)/2; #area of triangle

The cross product of the direction vectors formed from any two distinct pairs of the triangle's vertices is a normal vector (nv) to the plane (i.e., graph of the linear function) containing the triangle. Thus, the plane consists of all direction vectors emanating from any vertex of the triangle (I chose P) whose dot (or inner) product with nv is 0. The 2-Norm (or magnitude or length or Euclidean norm) of nv is twice the area of the triangle.

I think that the problem occurs because Maple is confused by possible situations where if certain parameters were equal then the geometric multiplicity of certain eigenvalues would change. A workaround for this is to change the parameters to "safe", "generic", and "independent" transcendental constants, compute the eigenvectors, then change the constants back to the parameters. Like this:

restart:

phi:= Pi:
A := <
          <0, 0, exp(I*k1) + m1, exp(I*k2) + m2>|
          <0, 0, exp(I*phi)*(exp(-I*k2) + m2), exp(-I*k1) + m1>|
          <exp(-I*k1) + m1, exp(-I*phi)*(m2 + exp(I*k2)), 0, 0>|
          <exp(-I*k2) + m2, exp(I*k1) + m1, 0, 0>
     >:
Subs:= [k1= exp(3), k2= exp(sqrt(2)), m1= sin(1), m2= cos(sqrt(2))]:
EV:= [LinearAlgebra:-Eigenvectors(eval(A, Subs))]:
eval(EV, (rhs=lhs)~(Subs)); simplify(%);

Doing this gives you two eigenvalues, each associated with two linearly independent eigenvectors. A quick visual scan suggests that these are the same as reported by Mathematica, although it may be a challenge to prove that conclusively. A starting point may be to have Mathematica convert its results to sin/cos form.

I substituted for all four parameters, but I suspect that it would still work with substituting for fewer.

First 113 114 115 116 117 118 119 Last Page 115 of 395