Carl Love

Carl Love

28110 Reputation

25 Badges

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

MaplePrimes Activity


These are answers submitted by Carl Love

For efficiency, it'd be better to use a custom procedure for this purpose rather than constructing a lengthy expression of ands and ors using seq. This is because in the vast majority of cases the evaluation can be terminated early without examining every (or even most) of the elements. Here's the procedure:

LexicographicOrder:= proc(a, b)
local k, na, nb;
   (na,nb):= (numelems(a), numelems(b));
   for k to min(na,nb) do
      if a[k] < b[k] then return true 
      elif a[k] > b[k] then return false
      end if
   end do;
   evalb(na < nb)
end proc:

If you're comparing lists of length 100, the procedure should be thousands of times faster for the average case than the construction with seqs. Constructing sequences of numerous subsequences is a relatively expensive operation in Maple.

Your second form and displayed form, both of which show the series as externally multiplied, are mathematical nonsense. But your first form, which has the series nested, makes sense. You said that it gives wrong answers. Note that the results below differ from your desired results by a factor of exactly 3. You must be missing a 3 somewhere.

restart:


SS:= sum(F[k-m]*sum(F[m-L]*sum(F[L-j]*F[j], j= 0..L), L= 0..m), m= 0..k):
eq:= (-1/(k+1))*(F[k] + sum((k-m+1)*F[k-m+1]*F[m], m= 0..k)/2+SS/20):
n:= 8:
F:= table():
for K from 0 to n do
   F[K+1]:= solve(eval(eq, k= K), F[K+1])
end do;

-(1/10)*F[0]^3-2

(3/200)*F[0]^2*(F[0]^3+20)

-(1/400)*F[0]*(F[0]^6+28*F[0]^3+160)

(7/16000)*F[0]^9+(63/4000)*F[0]^6+(3/20)*F[0]^3+1/5

-(9/800000)*F[0]^2*(7*F[0]^9+308*F[0]^6+4160*F[0]^3+16000)

(3/16000000)*F[0]*(77*F[0]^12+4004*F[0]^9+70080*F[0]^6+448000*F[0]^3+640000)

-(429/160000000)*F[0]^15-(1287/8000000)*F[0]^12-(3459/1000000)*F[0]^9-(777/25000)*F[0]^6-(171/1750)*F[0]^3-6/175

(9/17920000000)*F[0]^2*(1001*F[0]^15+68068*F[0]^12+1735776*F[0]^9+20137600*F[0]^6+100121600*F[0]^3+145920000)

-(1/179200000000)*F[0]*(17017*F[0]^18+1293292*F[0]^15+38152576*F[0]^12+542971520*F[0]^9+3722624000*F[0]^6+10304000000*F[0]^3+5836800000)

 

Download Solve_series.mw

Note that the -1/(k+1) coefficient in front of eq is superfluous. It couldn't possibly change the results because eq is equated to 0.

The polynomial equation is sixth degree in v. It's mathematically impossible to solve it. This is not a limitation of Maple.

The package linalg, and the commands evalm and crossprod had been superceded long before Maple 13. Also, you can reduce the repetition in your code by using elementwise operators:

restart:
ly:= 9.4607e15:
M0:= 1.99e30:

M:= M0*[2.20, 2.00, 1.50, 3.00]:

r:= [<-3, 3, 0>, <0, -2, 0>, <1, 2, 0>, <6, 4, 0>]*ly:
v:= [<25, 15, 0>, <20, -20, 0>, <-5, -25, 0>, <15, 0, 0>]*1e3:

Mtot:= `+`(M[]):

rcm:= `+`((M*~r)[])/Mtot:
vcm:= `+`((M*~v)[])/Mtot:

with(LinearAlgebra): 
add(r[i] &x (M[i]*v[i]), i= 1..nops(v)) - rcm &x (Mtot*vcm);

Other than the names, there is no similarity between eval and evalf. The command for evaluating an expression with free variables at values for those variables (a "point") is eval. The command for converting expressions with no free variables into complex decimal approximations is evalf. If you use evalf on an expression with free variables, then it is simply applied to the subparts that don't. It's not possible to use one command in place of the other.

You'll need to post the series problem that you're talking about to get further advice.

Induction is a proof technique. The corresponding programming technique is called recursion. This means that the procedure calls itself.

(A half a line of code omitted.)

That's all there is to it.

 

 

The error message says that it wants an indication of the dependent variables. So, change the pdsolve command to

pdsolve([pde, bc], u(x,t))

Diff(int(value(z), t), t);

The integrand has singularities at -I, 2*I, and 5. All of these are at distance greater than 1 from sqrt(11/2). So no singularities are enclosed by the path of integration. So the answer is 0.

If Maple ever gives you 2D output that you don't understand, just use lprint(%) to see the 1D equivalent. In this case, it will reveal the LerchPhi.

I don't think that dsolve's BVP solvers can handle complex numbers. But your BVP contains a parameter, A3. We can solve for this parameter and then use one of dsolve's IVP solvers, which do handle complex numbers. To solve for A3, note that the right side of your ODE doesn't contain the dependent variable, u. Thus we can use a combination of numeric integration and fsolve.
 

restart:

A1:= 5.5:  n:= .59:  A2:= 11818.:  h0:= 0.402e-3:
L:= .1:  dpx := -11823.9:  uc:= 0.44e-2:

ODE:= (A3,y)->
   (h0^(n+1)*L/sqrt(n)*(A1*exp(sqrt(n)*y/L)-A2*exp(-sqrt(n)*y/L))+dpx*y*h0+A3)^(1/n)
;

proc (A3, y) options operator, arrow; (h0^(n+1)*L*(A1*exp(sqrt(n)*y/L)-A2*exp(-sqrt(n)*y/L))/sqrt(n)+dpx*y*h0+A3)^(1/n) end proc

ODEINT:= proc(A3)
option remember;
local y;
   evalf(Int(ODE(A3,y), y= 0..1, epsilon= 1e-7)) - uc
end proc:

ReINT:= proc(A3x, A3y)
   Digits:= 15:
   Re(ODEINT(A3x + I*A3y))
end proc:

ImINT:= subs(Re= Im, eval(ReINT)):

Digits:= 7:
a3:= fsolve([ReINT, ImINT]);

[2.3776470627663, -1.0154578732062]

A3:= Complex(a3[]);

2.3776470627663-1.0154578732062*I

Solve as IVP:

Digits:= 15:
sol:= dsolve({diff(u(y),y) = ODE(A3,y), u(0)=0}, numeric):

plots:-odeplot(
   sol, [[y, Re(u(y))], [y, Im(u(y))]], y= 0..1,
   legend= [real, imag], labels= [y, u(y)]
);

Verify that boundary condition at u(1) is satisfied:

abs(eval(u(y), sol(1)) - uc);

0.489191321784278e-7

sol(.5);

[y = .5, u(y) = .504669074500681-.965536950733450*I]

``

Download Complex_BVP.mw

 

The cause of the anomaly that you described is that print doesn't return values; it merely displays them in your worksheet. The purpose of print is to display supplementary information from the middle of a procedure. It is not appropriate to use it at the end of a procedure to return the procedure's main return value. If you do use it like this, the actual returned value (not the displayed value) is NULL, which is of type exprseq, which is what you're experiencing. So, all you need to do is change print(Hm) to simply Hm.

I think that you may have the mistaken idea that a "statement" needs to be on a single line. Actually, there's no practical limit to the number of lines that it can take, the number of substatements internal to it, or the depth of their nesting.

Change

if m1 <= 2^(2*k-1) then t1 := (C+(-A1*(m1^2)))/A2 end if;
Do (t1);

to

if m1 <= 2^(2*k-1) then
   t1 := (C+(-A1*(m1^2)))/A2;
   Do (t1)
end if;

And do likewise for the other if-then blocks.

The line

E2:= V(x2);

should be

E2:= V(x2(i));

 

It's not a source of the error, but I'm concerned about your use of Seed. What is its significance? Why not use rand()?

You should not use randomize() until you have your code running correctly without it because its use will make some bugs harder to find.

Please post code in plaintext form, or post a worksheet. I was unable to test my answer above because I couldn't run your code because it's in a form that can't be copied-and-pasted.

If the function has a Taylor series at x=a, then that's what's returned. For complex-valued functions (which most functions are), having a Taylor series at x=a is equivalent to the function being differentiable throughout an open set containing x=a. Next, if it has a Laurent series at x=a, then that's what's returned. Having a Laurent series at x=a is equivalent to there being a positive integer n such that (x-a)^n*f(x) has a Taylor series at x=a. Finally, if a more-generalized type of series can be found, then that's what's returned. I don't know the details of these more-generalized series, but my first guess is that they're always of the form M((x-a)*g(x-a)), where M(z) is a Maclaurin series and g(x) is a relatively simple function (such as ln(x) or x^r for some r, 0 < |r| < 1) that is not differentiable at x=a.

First 193 194 195 196 197 198 199 Last Page 195 of 395