Carl Love

Carl Love

28055 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

I can't look at your worksheet because MaplePrimes doesn't accept the .maple format. If you resave and repost it as a .mw, I'll be able to look at it. This should be no problem because your work shouldn't require any auxilliary data.

There are two errors in the printed statement of the problem:

  • and Re are used to represent the same thing. I just changed Re to R.
  • The first initial condition for the third equation should be F2(0)=0, not 
    F0(0)=0.

Here's my solution, which matches what's stated in the paper:


 

Eq0:= {diff(F0(eta), eta$3) = 0, F0(0)=1, D(F0)(0)=0, (D@@2)(F0)(0)=a}:

Sol0:= dsolve(Eq1);

F0(eta) = (1/2)*a*eta^2+1

(1)

Eq1:= {
   diff(F1(eta), eta$3) +
      2*alpha*R*F0(eta)*diff(F0(eta),eta) +
      alpha^2*(4-H)*diff(F0(eta),eta)
   = 0,
   F1(0)=0, D(F1)(0)=0, (D@@2)(F1)(0)=0
}:

Sol1:= dsolve(eval(Eq1, Sol0));

F1(eta) = a*alpha*(-(1/120)*R*a*eta^6+(1/24)*(H*alpha-2*R-4*alpha)*eta^4)

(2)

Sol1:= collect(Sol1, eta);

F1(eta) = -(1/120)*a^2*alpha*R*eta^6+a*alpha*((1/24)*H*alpha-(1/12)*R-(1/6)*alpha)*eta^4

(3)

Eq2:= {
   diff(F2(eta), eta$3) + alpha^2*(4-H)*diff(F1(eta),eta) +
      2*alpha*R*F0(eta)*diff(F1(eta),eta) +
      2*alpha*R*F1(eta)*diff(F0(eta),eta)
   = 0,
   F2(0)=0, D(F2)(0)=0, (D@@2)(F2)(0)=0
}:      

Sol2:= dsolve(eval(Eq2, {Sol0, Sol1}));

F2(eta) = (1/30)*a*alpha^2*((1/360)*R^2*a^2*eta^10+(1/336)*(-9*H*R*a*alpha+18*R^2*a+36*R*a*alpha)*eta^8+(1/120)*(5*H^2*alpha^2-20*H*R*alpha-40*H*alpha^2+20*R^2+80*R*alpha+80*alpha^2)*eta^6)

(4)

Sol2:= collect(Sol2, eta);

F2(eta) = (1/10800)*a^3*alpha^2*R^2*eta^10+(1/30)*a*alpha^2*(-(3/112)*H*R*a*alpha+(3/56)*R^2*a+(3/28)*R*a*alpha)*eta^8+(1/30)*a*alpha^2*((1/24)*alpha^2*H^2-(1/6)*alpha*H*R-(1/3)*alpha^2*H+(1/6)*R^2+(2/3)*R*alpha+(2/3)*alpha^2)*eta^6

(5)

 


 

Download SequentialDsolve.mw

 

The command is currentdir. Example:

currentdir("C:/Users/carl/desktop");

The forward slash / can always be used as the directory separator in Maple filenames, even on Windows.

A distinction needs to be made between equations and assignments. This distinction can be a slightly difficult concept to grasp for new programmers. It is often the first nontrivial concept one encounters when learning a first computer language. This distinction exists in all computer languages, but it is more pronounced in symbolic languages such as Maple, where variables can be meaningful even when they don't have values (we say that such variables are symbolic).

Equations are declarative statements (in the logical/linguistic sense). They can be true or false depending on the values of their variables. With symbolic variables, "truth" can be more subtle: An equation may be true for all, for some, or for no values of its variables. In Maple, equations are expressed with the equals sign =.

On the other hand, assignments are imperative statements (i.e., commands) that force a variable to be a certain value. (That certain value may possibly still contain symbolic variables.) Thus, there is no question of their truth. In Maple, assignments are expressed with := or occasionally with the assign or unassign commands. Assignments are often expressed in mathematical writing with "let", as in "Let x = 3."

Your problem would be best done using equations. Like Acer, I disdain making assignments to the directly meaningful variables of a problem---those variables that are meaningful outside the context of a computer program. 

Here is the solution to your problem. Note how conveniently the Solution is expressed: It shows the numeric values of all variables in one place. These numbers are not assigned to the variables. Rather, solve is saying that if the variables had those values, the equations would be true.

System:= {
   totalsalesx = 60.1 + tax + profit,
   tax = 0.3*profit,
   profit = 0.1*totalsalesx
}:
Solution:= solve(System);
           {profit = 6.908045977, tax = 2.072413793, totalsalesx = 69.08045977}

 

Here is a procedure for it:

DroppedContourPlot3d:= proc(
   Z::name= algebraic, #function to plot
   X::name= range(realcons), #x-axis-parameter range
   Y::name= range(realcons), #y-axis-parameter range
   {
      #fraction of empty space above and below plot (larger "below"
      #value improves view of dropped-shadow contourplot):
      zmargin::[realcons,realcons]:= [.05,0.15],
      contouropts::list({name, name= anything}):= [],
      surfaceopts::list({name, name= anything}):= []    
   }
)
option `Author: Carl Love <carl.j.love@gmail.com> 19-Nov-2019`;
local
   LX:= lhs(X), RX:= rhs(X), LY:= lhs(Y), RY:= rhs(Y),
   C:= plots:-contourplot(
      rhs(Z), X, Y,
      #These default plot options can be changed by setting 'contouropts': 
      'contours'= 5, 
      'filled', 'coloring'= ['yellow', 'orange'], 'color'= 'green', 
      contouropts[]
   ),
   P:= plot3d(
      rhs(Z), X, Y,
      #These default plot options can be changed by setting 'surfaceopts': 
      'style'= 'surfacecontour', 'contours'= 6,
      surfaceopts[]
   ),
   U, L #z-axis endpoints after margin adjustment
;
   #Stretch z-axis to include margins:
   (U,L):= ((Um,Lm,M,m)-> (M*(Lm-1)+m*Um, M*Lm+m*(Um-1)) /~ (Um+Lm-1))(
      zmargin[],
      (max,min)(op(3, indets(P, 'specfunc'('GRID'))[])) #actual z-axis range
   );
   plots:-display(  #final plot to display/return
      [
         plots:-spacecurve(
            { 
               [  #(y,z)-plane backwall:
                  [lhs(RX),rhs(RY),U], [rhs(RX),rhs(RY),U],
                  [rhs(RX),rhs(RY),L]
               ],
               [  #((x,z)-plane backwall:
                  [rhs(RX),rhs(RY),U], [rhs(RX),lhs(RY),U],
                  [rhs(RX),lhs(RY),L]
               ] 
            },
            'color'= 'grey', 'thickness'= 0
         ),
         plottools:-transform((x,y)-> [x,y,L])(C), #dropped contours
         P
      ],
      #These default plot options can be changed simply by putting the option
      #in the DroppedCountourPlot3d call:
      'view'= ['DEFAULT', 'DEFAULT', L..U], 'orientation'= [-135, 75], 
      'axes'= 'frame',
      'labels'= [lhs(X), lhs(Y), 'z'], 
      'labelfont'= ['TIMES', 'BOLDOBLIQUE', 16],
      'axesfont'= ['TIMES', 12],
      'projection'= 2/3,   
      _rest
   )
end proc
: 

Example usage:

DroppedContourPlot3d(z=x^2+y^2, x= -2..2, y= -2..2);

I didn't plot the same function that you specified because it's too difficult to read your unformatted Question, but it should be straightforward to plot it as above.

My procedure automatically chooses a plane for the 2d plot to optimize the view; it doesn't necessarily use plane z=0.

Since you tried partial fractions, I thought that maybe you'd want a step-by-step work through of the "telescoping" method by which the sum of such a series is found. Here it is.
 

restart:

T:= 1/i/(i+1)/(i+2)/(i+3);

1/(i*(i+1)*(i+2)*(i+3))

(1)

We want this sum:

S:= Sum(T, i= 1..infinity);

Sum(1/(i*(i+1)*(i+2)*(i+3)), i = 1 .. infinity)

(2)

We suspect "telescoping" when the general term contains parts that also occur in subsequent terms, i.e., i+1, i+2, i+3. An analysis of it usually begins by breaking the general term itself into a sum. In this case, that means a partial-fraction decomposition.

F:= unapply(convert(T, parfrac, i), i);

proc (i) options operator, arrow; (1/6)/i+(1/2)/(i+2)-(1/6)/(i+3)-(1/2)/(i+1) end proc

(3)

Add an arbitrary number of consecutive terms from anywhere in the series to see the effect of "telescoping" cancellations.

sum(F(i+j), j= 0..n);

-(1/6)/(i+n+1)+(1/3)/(i+n+2)-(1/6)/(i+n+3)+(1/6)/i-(1/3)/(i+1)+(1/6)/(i+2)

(4)

This is called a partial sum of the series. Its limit is the sum of the series. Obviously the terms with n in the denominator go to 0.

limit(%, n= infinity);

(1/3)/(i*(i+1)*(i+2))

(5)

The actual sum that we want begins at i=1.

eval(%, i= 1);

1/18

(6)

Check if Maple concurs:

value(S);

1/18

(7)

 


 

Download TelescopingSeries.mw

Hints:

1. The command fsolve when applied to a univariate polynomial and given the option complex will return all the roots. 

2. The command abs returns the modulus.

3. The command Digits:= 15 should give you enough leeway that you can fully trust the first 10 digits.

In writing the assignment, your instructor inadvertently used f in two different ways: as a coefficient and as a function. This is of course confusing to Maple. Tom's Answer gets around this by using and h as the functions. 

One way to interpolate a data list for a smoother plot is

S:= [seq(j^2, j= -5..5)]:
plot(Interpolation:-SplineInterpolation([$1..nops(S)], S), 1..nops(S));

The SplineInterpolation command converts the data from a list to a function, thus it can be plotted with plot but not listplot.

You mentioned "time series". If you are working with data that are actually indexed by absolute times (as measured by calendars and clocks), then look at Maple's TimeSeries package.
 

Here are the first and last steps to do it. There's a step in the middle that I can't be definite about without seeing a worksheet containing the tensor. You can post your worksheet by editing your Question and using the green uparrow on the toolar. You edit the Question by using the "More..." pull-down menu in its lower right corner.

First step: Create an Array of indexed names

Arr:= Array((1..4)$6, ()-> index(myT, args)):

This step creates an array of 4096 names with indices in 6 dimensions, each running from 1 to 4. The names are all of the form myT[i,j,k,l,m,n].

Middle step: Match the names to the tensor's components: I'll suppose that the tensor is named T. This step may be as easy as

Matches:= Arr =~ T: ,

but I'm not sure if the elementwise operator =~ will work in your case. That's why I need to see your worksheet.

Last step: Convert to a set

mySet:= {seq(Matches)}:

Now mySet is set of 4096 equations, each of the form 
myT[i,j,k,l,m,n] = ....

Example: Retrieving a value: I don't recommend that you actually assign the values (although you may do it with assign~(mySet)). Rather, any desired value can be retrieved as in this example:

eval(myT[1,2,3,4,1,2], mySet);

plot(ListTools:-Enumerate(S));

The shortest command I can think of for this is

plot(<<[$1..nops(S)]>|<S>>);

If you're willing to forego any special handling of exceptional conditions (such as fsolve not being able to find a solution), and you're satisfied with just one solution, then it can be done as a single command:

Z:= eval(Matrix((2,2), symbol= z), fsolve({seq(A)}, indets(A,name)=~ 0..1)); 

Your dsolve command is missing a closing brace }. This is the type of easy-to-find error to look for when you get a "unable to match delimiters" error message.

My previous Answer only needs a trivial modification for the new integrals. You seem to have completely ignored that Answer. Here is the modification:

U:= Matrix(
   N, 
   (i,j)-> 
      int(
         psi[i](x)*omega[iquo(i-1,M)+1](x)*
            int(u(x,t)*psi[j](t)*omega[iquo(j-1,M)+1](t), t= 0..1),
         x= 0..1
      )
);

 

It's fairly easy to make them wider, as ellipses, by padding the labels before and after with spaces:

G:= GraphTheory:-CompleteGraph(5):
G1:= GraphTheory:-RelabelVertices(
   G, 
   sprintf~("    %a    ", GraphTheory:-Vertices(G))
):
GraphTheory:-DrawGraph(G1, stylesheet= "legacy");

Thanks. Your Matlab work shows both sufficient effort and sufficient understanding that I'm willing to share my Maple code.

Your solution of the differential equation is correct. The plotting "blip" that you showed (I wish that you hadn't deleted it!!) comes, I think, from Matlab's spline command. I don't have Matlab to check.
 

restart:

Digits:= 15:

#
#All problem specs are in this section.
#
params:= {X0= 0.05, B= 3.9, Nrange= 0..30, Yrange= 0..1}:
params:= params union
   #recurrence formula:
   {R= subs(params, n-> B*R(n)*(1-R(n)))}
:
(Nrange, Yrange):= (0..30, 0..1):

#
#Everything else is generic (and could be applied to any 1st-order recurrence
#once the parameters are specified).
#
 
#recursive procedure for the difference equation:
X:= subs[eval](
   _R= eval(R, params), R= thisproc,
   proc(n::nonnegint) option remember; _R(n-1) end proc
):
X(lhs(Nrange)):= eval(X0, params) #Put initial value in remember table.
:

#equivalent differential equation (which has a simple symbolic solution in
#this case):
F:= rhs(
   dsolve(
      subs(R= Y,  eval({diff(Y(n),n)+Y(n) = R(n), Y(lhs(Nrange))=X0}, params))
   )
);

29/(39+541*exp(-(29/10)*n))

#plot of both
plot(
   [F, [seq([k,X(k)], k= Nrange)]], n= Nrange,
   style= [line, pointline], thickness= [1,0], size= [600,400],
   axes= boxed, view= [Nrange, Yrange],
   legend= [continuous, discrete], gridlines,
   title= typeset(
      "Logistic example:  ", beta=eval(B,params), ", ", 'X'[0]=eval(X0,params)
   ),
   titlefont= [TIMES, BOLD, 16],
   labels= [`n\n`, 'X'[n]], labelfont= [HELVETICA, BOLD, 14]
);

 


 

Download LogisticMap.mw

Please let me know if you have any trouble with the above code, or any other questions. I don't have Maple 18 to check this, but I didn't conciously use anything that wouldn't work in Maple 18. Note that the plot looks much better in the actual worksheet than it's rendered above by MaplePrimes.

First 125 126 127 128 129 130 131 Last Page 127 of 395