Carl Love

Carl Love

28050 Reputation

25 Badges

12 years, 336 days
Himself
Wayland, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity


These are replies submitted by Carl Love

@abiodunowoyemi9 It's a version issue. You should use the pull-down list in the Question header to state which Maple version you are using. Below, I've removed some of the newer features from my code. Please try this:
 

Direct application of Method of Undetermined Coefficients

restart
:

V:= [v__||(0..5)]:  v||(0..5):= V[]:  Vs:= ()
:

#right sides of the 6 odes:
RHS:= -[
    0, v0^3, 3*v0^2*v1, 3*v0*v1^2 + 3*v0^2*v2,
    6*v0*v1*v2 + 3*v0^2*v3 + v1^3,
    6*v0*v1*v3 + 3*v1^2*v2 + 3*v0^2*v4 + 3*v0*v2^2
]:
print~(RHS)
:

0

-v__0^3

-3*v__0^2*v__1

-3*v__0^2*v__2-3*v__0*v__1^2

-3*v__0^2*v__3-6*v__0*v__1*v__2-v__1^3

-3*v__0^2*v__4-6*v__0*v__1*v__3-3*v__0*v__2^2-3*v__1^2*v__2

ICs:= [[1,0], [0,0]$5] #initial conditions of the 6 odes
:

for k to nops(RHS) do
    S:= C1*sin(t) + C2*cos(t); #general homogenous solution of v'' + v = 0
    #Perform method of undetermined coefficients termwise:
    rs:= combine(eval(RHS[k], [Vs]));
    for f in indets(rs, specfunc({sin,cos})) do #for each term on right side
        C:= coeff(rs, f);
        d:= degree(C,t);
        o:= op(f);
        P:= add(a[i]*t^i, i= 0..d)*sin(o) + add(b[i]*t^i, i= 0..d)*cos(o);
        if o=t then P:= P*t fi; #extra power of t   
        S:= S + eval(P, solve(identity(diff(P,t$2)+P=C*f, t)))
    od;
    #Apply initial conditions to determine C1 and C2:
    Vs:= Vs, V[k] = (combine@eval)(S, solve(eval([S, diff(S,t)], t=0)=~ICs[k]))
od
:
print~([Vs])
:

v__0 = cos(t)

v__1 = -(1/32)*cos(t)-(3/8)*sin(t)*t+(1/32)*cos(3*t)

v__2 = -(9/128)*cos(t)*t^2+(3/32)*sin(t)*t-(9/256)*t*sin(3*t)+(1/1024)*cos(5*t)-(3/128)*cos(3*t)+(23/1024)*cos(t)

v__3 = (9/1024)*sin(t)*t^3-(81/4096)*t^2*cos(3*t)+(135/4096)*cos(t)*t^2-(15/8192)*t*sin(5*t)+(279/8192)*t*sin(3*t)-(207/4096)*sin(t)*t+(1/32768)*cos(7*t)-(3/2048)*cos(5*t)+(297/16384)*cos(3*t)-(547/32768)*cos(t)

v__4 = (27/32768)*cos(t)*t^4-(99/16384)*sin(t)*t^3+(243/32768)*t^3*sin(3*t)-(1359/65536)*cos(t)*t^2+(1539/65536)*t^2*cos(3*t)-(225/131072)*t^2*cos(5*t)-(21/262144)*t*sin(7*t)+(8997/262144)*sin(t)*t-(3915/131072)*t*sin(3*t)+(825/262144)*t*sin(5*t)-(9/131072)*cos(7*t)+(6713/524288)*cos(t)-(15121/1048576)*cos(3*t)+(883/524288)*cos(5*t)+(1/1048576)*cos(9*t)

v__5 = -(81/1310720)*sin(t)*t^5-(783/1048576)*cos(t)*t^4+(2187/1048576)*t^4*cos(3*t)+(1125/1048576)*t^3*sin(5*t)+(4635/1048576)*sin(t)*t^3-(10935/1048576)*t^3*sin(3*t)+(6975/2097152)*t^2*cos(5*t)+(15777/1048576)*cos(t)*t^2-(441/4194304)*t^2*cos(7*t)-(96795/4194304)*t^2*cos(3*t)-(16575/4194304)*t*sin(5*t)-(108357/4194304)*sin(t)*t+(1659/8388608)*t*sin(7*t)+(108243/4194304)*t*sin(3*t)-(27/8388608)*t*sin(9*t)-(921/524288)*cos(5*t)-(42397/4194304)*cos(t)+(1757/16777216)*cos(7*t)-(3/1048576)*cos(9*t)+(394701/33554432)*cos(3*t)+(1/33554432)*cos(11*t)

plot(
    rhs~([Vs]), t= 0..15,
    color= [red, green, blue, yellow, cyan, black],
    axes=boxed, gridlines=false, legend= V
);

#Compare with numeric solution:
sys:= [
    (diff~(V(t), t$2) +~ V(t) =~ RHS(t))[],
    (V(0)=~ICs[..,1])[], (D~(V)(0)=~ICs[..,2])[]
]:

sol:= dsolve(sys, numeric):

plots:-odeplot(
    sol, `[]`~(t, V(t)), t= 0..15,
    color= [red, green, blue, yellow, cyan, black],
    axes=boxed, gridlines=false, legend= V
);

 

Download UndeterminedCoeffs.mw

@vv Yes, I noted that and replied to him above.

@tomleslie You have exactly the same typos on the right sides of the 3rd, 4th, and 5th equations as VV has. I guess he used your encoding of the equations.

@vv You have a few typos on the right-sides of the 3rd, 4th, and 5th equations. I've corrected these below. Running the corrected code (in Maple 2019), I experience the same thing that I said above: an undetermined very long time for the 6th solution.

restart;
sys:= [ 
[diff(v__0(t),t$2)+v__0(t)=0,  v__0(0)=1, D(v__0)(0)=0],
[diff(v__1(t),t$2)+v__1(t)=-v__0(t)^3,  v__1(0)=0,  D(v__1)(0)=0],
[diff(v__2(t),t$2)+v__2(t)=-3*v__0(t)^2*v__1(t),  v__2(0)=0,  D(v__2)(0)=0],
[diff(v__3(t),t$2)+v__3(t)=-3*v__0(t)*v__1(t)^2-3*v__0(t)^2*v__2(t),  v__3(0)=0, D(v__3)(0)=0],
[diff(v__4(t),t$2)+v__4(t)=-6*v__0(t)*v__1(t)*v__2(t)-3*v__0(t)^2*v__3(t)-v__1(t)^3,  v__4(0)=0,  D(v__4)(0)=0],
[diff(v__5(t),t$2)+v__5(t)=-6*v__0(t)*v__1(t)*v__3(t)-3*v__1(t)^2*v__2(t)-3*v__0(t)^2*v__4(t)-3*v__0(t)*v__2(t)^2,  v__5(0)=0, D(v__5)(0)=0]
]:
print~(sys);
s:=NULL:
for i to 6 do  s:=s, dsolve(eval(sys[i], [s])); print([s][-1]) od:

 

@tomleslie It was not my intention to restrict n to integers. Yes, if that had been my intention, I would've included adaptive= false. My intention was to make sure that every integer was included along with whatever other points plot decided to include on its own.

Yes, I understand your point about the inner floors. They aren't needed in my Answer. I mistakenly thought that I was using the OP's original function, but I was actually using your modification of it due to a hasty copy and paste. Sorry about that. So, the code of my Answer should be changed to 

f:= n-> ceil(sqrt(4*n)) - floor(sqrt(2*n)) - 1:
plot(f, 10...100, sample= [$10..100]);

The distance function that I gave before should be modified to account for cases that require no horizontal travel. This doesn't affect the final results previously obtained because your example doesn't have any such pairs of points. The new distance function is 

dist:= proc(p::[integer,integer], q::[integer,integer])
local H;
    if p[1]=q[1] then abs(p[2]-q[2])
    else min((add@abs~)~([seq]([p[2]-H, p[1]-q[1], H-q[2]], H= [1,10])))
    fi
end proc      
:

 

@jalal In the line defining fnew, try changing ^+ to ^%T. They mean the same thing---transpose the matrix---but ^+ only works in (plaintext) 1D input.

@Teep Here is the code to highlight with arrows the direction of travel for the solution to the traveling salesman problem.

restart:
GT:= GraphTheory: Pl:= plots: Co:= combinat:
dist:= proc(p::[integer,integer], q::[integer,integer])
local H;
    min((add@abs~)~([seq]([p[2]-H, p[1]-q[1], H-q[2]], H= [1,10])))
end proc      
:
HighlightTour:= proc(    
    #edges:
    E::set([And(set({name, string, integer}), 2 &under nops), numeric]),
    #tour:
    T::list({name, string, integer}), #the  tour
    #vertices (optional, for positions only):
    V::list({name, string, integer} = [realcons,realcons]):= 0
)
uses GT= GraphTheory;
local G0,Gt;
    (G0,Gt):= GT:-Graph~(
        [selectremove](e-> e[1]::set, subs((`{}`=`[]`)~(T[..-2],T[2..]), E)),
        `if`(V=0, [][], 'vertexpositions'= rhs~(V))
    )[];
    GT:-HighlightTrail(Gt, T, _rest, 'inplace');
    G0,Gt
end proc
:
pts2 := [A = [5, 0], B = [8, 9], C = [2, 8], D = [4, 3], E = [7, 4]]:
Edgs:= {seq}([({op}@lhs~)(e), dist(rhs~(e)[])], e= Co:-choose(pts2, 2)):
G1:= GT:-Graph(Edgs):
(distance, Tour):= GT:-TravelingSalesman(G1);
            distance, Tour := 38, [A, D, C, B, E, A]
(G0,Gt):= HighlightTour(Edgs, Tour, pts2);
G0, Gt := Graph 1: an undirected weighted graph with 5 vertices 
   and 5 edge(s), 
  Graph 2: a directed weighted graph with 5 vertices and 5 arc(s)

Pl:-display(GT:-DrawGraph~([G0,Gt], style= circle, stylesheet= "legacy"));

#Using actual vertex positions:
Pl:-display(GT:-DrawGraph~([G0,Gt], stylesheet= "legacy"));

To see the tour (red) edges only, remove G0 from either of those plot commands.
 

I'm having some trouble at the moment loading the workbook format (*.maple). Would it be possible for you to duplicate the error in the worsheet (*.mw) format?

If you don't give any constraints, then obviously there's no upper bound. For example, you could make A1 any large value and make all the other variables 0.

@Vini Your new Question, which I just deleted twice, is too close to this Question. Please discuss the matter in this thread.

At this point, I don't see any mathematical difference whatsoever in the two Questions; so, please try to say more precisely what you think the difference is.   

@Carl Love The above can be improved a little. I don't like that the 17 bars are crammed into interval 0..16. The interval should be -0.5..16.5. The improvement is a command called tallyinto.

B:= 16:
plots[display]([
    stats[statplots, histogram](
        stats[transform, tallyinto](
            map(`+`@op@convert, [$0..2^B-1], base, 2),
            [seq](k-1/2..k+1/2, k= 0..B)
        ),
        area= 2^B, color= pink, thickness= 0
    ),        
    plot(binomial(B,x), x= 0..B, color= blue, thickness= 2)
]);

@logan35 Oh my, Maple 7. That's a shame. In general, it's going to be difficult for you to get usable answers here on MaplePrimes. Very few of us have versions that old (about 20 years). Luckily, I've been using Maple long enough to have some memory of the old commands. So, try this:

plots[display]([
    stats[statplots, histogram](
        map(`+`@op@convert, [$0..2^16-1], base, 2), 
        area= 2^16, numbars= 17, color= pink
    ),        
    plot(binomial(16,x), x= 0..16, color= blue, thickness= 3)
]);

 

@logan35 The entire program is simply the one line that I gave:

Statistics:-Histogram((`+`@op@Bits:-Split)~([$0..2^16]), binwidth= 1);

@Teep Sorry, that was my fault. Try again; I just edited the code.

First 92 93 94 95 96 97 98 Last Page 94 of 709