Carl Love

Carl Love

28065 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

If you have a working plot3d command such as 

plo3d(u, x= -3..3, t= -3..3);

then simply change it to

plots:-contourplot(u, x= -3..3, t= -3..3);

Note that I use u rather than [u]. The [u] is acceptable to plot3d but not to contourplot.

Here is the solution with plots for the equation solved as a BVP and solved via HPM (Homotopy Perturbation Method). I didn't make the table, but you could extract the data for the table from the matrices in the plots (indeed, for 201 values of eta in 0..1) just like I did in the code below. To my mind, all that matters is the maximum error over the interval 0..1, so that's what I computed. I think that the table was just filler for a relatively short paper.
 

restart #Remove restart if using Maple 2019.2 due to a serious bug.
:

All Equation Eq numbers refer to the corresponding Equations in the paper.

 

The fundamental BVP:

Eq6:= diff(f(eta),eta$3) +
    alpha*(2*R*f(eta) + (4-H)*alpha)*diff(f(eta),eta)
= 0;
Eq7:= f(0)=1, D(f)(0)=0, f(1)=0:
BVP:= {Eq6, Eq7}:

#Parameters:
params:= {alpha= 7.5*Pi/180, R= 50}:
#Note that I converted alpha to radians!
#R is "Re" in the paper.

Hlist:= [0, 250, 500, 1000]: #Hartmann numbers

diff(diff(diff(f(eta), eta), eta), eta)+alpha*(2*R*f(eta)+(4-H)*alpha)*(diff(f(eta), eta)) = 0

for k to nops(Hlist) do
    Ha:= Hlist[k];
    Sol_f[Ha]:= dsolve(eval(BVP, params union {H= Ha}), numeric);
    P[Ha]:= plots:-odeplot(
       Sol_f[Ha], [eta, f(eta)], legend= [H=Ha],
       color= COLOR(HUE, (k-1)/(nops(Hlist)-1)*.85)
    )
od:

#Fig. 2:
plots:-display(
   [seq(P[k], k= Hlist)], gridlines= false,
   title= "Velocity profiles, solved as BVPs"
);

Homotopy Perturbation Method (HPM)

Eq17:= {diff(F0(eta), eta$3) = 0, F0(0)=1, D(F0)(0)=0, (D@@2)(F0)(0)=a};
Eq18:= {
   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
};
Eq19:= {
   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
};

{F0(0) = 1, diff(diff(diff(F0(eta), eta), eta), eta) = 0, (D(F0))(0) = 0, ((D@@2)(F0))(0) = a}

{diff(diff(diff(F1(eta), eta), eta), eta)+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}

{diff(diff(diff(F2(eta), eta), eta), eta)+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}

Eq20:= dsolve(Eq17);

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

Sol:= dsolve(eval(Eq18, Eq20)):
Eq21:= collect(Sol, eta);

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

Sol:= dsolve(eval(Eq19, {Eq20, Eq21})):
Eq22:= collect(Sol, eta);

F2(eta) = (1/10800)*alpha^2*a^3*R^2*eta^10+(1/30)*alpha^2*a*(-(3/112)*alpha*H*R*a+(3/56)*R^2*a+(3/28)*R*alpha*a)*eta^8+(1/30)*alpha^2*a*((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

f_HPM:= unapply(eval(F0(eta)+F1(eta)+F2(eta), {Eq20,Eq21,Eq22}), eta);

proc (eta) options operator, arrow; (1/2)*a*eta^2+1-(1/120)*alpha*a^2*R*eta^6+alpha*a*((1/24)*H*alpha-(1/12)*R-(1/6)*alpha)*eta^4+(1/10800)*alpha^2*a^3*R^2*eta^10+(1/30)*alpha^2*a*(-(3/112)*alpha*H*R*a+(3/56)*R^2*a+(3/28)*R*alpha*a)*eta^8+(1/30)*alpha^2*a*((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 end proc

Solve for a:

a_sol:= solve({f_HPM(1)=0}, {a})[1]: #[1] is only real solution.

for k to nops(Hlist) do
    Ha:= Hlist[k];
    PH[Ha]:= plot(
       eval(eval(f_HPM(eta), a_sol), params union {H= Ha}),
       eta= 0..1, legend= [H=Ha],
       sample= [seq(op([1,1], P[Ha])[..,1])], adaptive= false,
       color= COLOR(HUE, (k-1)/(nops(Hlist)-1)*.85)
    )
od:

#Fig. 2:
plots:-display(
    [seq(PH[k], k= Hlist)], gridlines= false,
    title= "Velocity profiles, solved via HPM",
    labels= ['eta', 'f']
);

Compute maximum difference between the two solutions for each H.

for Ha in Hlist do
    err[` H`=Ha]:= LinearAlgebra:-Norm(
        Vector(op([1,1], P[Ha])[..,2] - op([1,1], PH[Ha])[..,2]),
        infinity #max norm
    )
od;

.129520180880002722

0.240134111942442163e-1

0.149072751558998462e-2

0.698175154954450150e-2

 


 

Download BVPvsHPM.mw

If your goal is simply to represent objects defined in real vector spaces with more than 3 dimensions, you can use animation and/or color. Animation makes time a dimension. There are several schemes (RGB, HSV, etc.) for squeezing in 3 dimensions of color. 

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.

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