sand15

1648 Reputation

16 Badges

11 years, 208 days

MaplePrimes Activity


These are replies submitted by sand15

@vv 

(y^2+1)/(x^4+y^4+1)

should be

(y^2+1)/(x^4+y^4+2)

After correction your result is  1.53978589107117

@Rouben Rostamian
vv accounted for the 'R' term in 

eval(f, [x=R*cos(t), y=R*sin(t)])*R

@nm 

Note that in the excerpt you provide it is written "We obviously have... u(∞)=0. By the latter we clearly mean lim(u(x), x=+∞)=0." ... which is nothing more than the author's confession of using a notational convention to achieve a more compact notation that is not mathematically rigorous.

You write later "I am here talking about using a CAS to solve an ode. Not about how one would write things for a math paper".
This is unbelievable: I get the impression that you would like Maple to solve a math problem that you haven't formulated correctly from a mathematical standpoint, all that just to ask you next why the solution isn't correct.
If you don't want to bother with mathematical rigor, you should try using some kind of AI to ask it how to solve your poorly formulated problem using Maple: perhaps it will interpret correctly some abusive notations and give you the way to proceed.

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

# You wrote

ode:=diff(y(x),x$2)-y(x)=0;

IC:= y(0)=1,limit(y(x),x=infinity)=0;

sol:=dsolve([ode,IC])

 

diff(diff(y(x), x), x)-y(x) = 0

 

y(0) = 1, limit(y(x), x = infinity) = 0

 

Error, (in dsolve) invalid input: op expects 1 or 2 arguments, but received 0

 

# I write

ode:=diff(y(x),x$2)-y(x)=0;

IC:= y(0)=1, y(L) = 0;

sol:=dsolve([ode,IC]);
limit(sol, L = +infinity)

diff(diff(y(x), x), x)-y(x) = 0

 

y(0) = 1, y(L) = 0

 

y(x) = exp(L)*exp(-x)/(-exp(-L)+exp(L))-exp(-L)*exp(x)/(-exp(-L)+exp(L))

 

y(x) = exp(-x)

(2)

# Same limits whatever the path taken to {x -> +infinity, L -> + infinity, x < L}

limit(limit(eval(sol, x = t*L), t = 1, left), L = +infinity);
limit(limit(eval(sol, x = t*L), L = +infinity), t = 1, left);

limit(limit(eval(sol, x = t*L), t = a, left), L = +infinity);
eval(%) assuming a > 0, a < 1

limit(y(L), L = infinity) = 0

 

limit(y(L), L = infinity) = 0

 

limit(y(a*L), L = infinity) = limit(exp(L)*exp(-a*L)/(-exp(-L)+exp(L))-exp(-L)*exp(a*L)/(-exp(-L)+exp(L)), L = infinity)

 

limit(y(a*L), L = infinity) = 0

(3)
 

 

Download My_point_of_view.mw

Note: had you interested in numerical solution of odes (pdes) with boundary conditions at infinity, you should chose some right boundary L >> C (usually based on a dimension analysis of these odes/pdes and the identification of some characteristic length C) and solve these odes/pdes on a bounded domain with extension L in the "infinite" direction.
Next, to ensure your choice of L wascorrect, you should solve again this same poblem for a larger value of L and verify its solutions is almost the same of the one you got earlier.
Another way, more technical bur lore correct, is to reformulate your problem using the Integral Equations Method (largely used in electromagnetism for instance).

What I did above, using the lright boundary L and making it tend to infinity, is nothing but what any people solving numerically this ode would do (solve it for L = L1 and L = L2 > L1, look at the differences between the solutions and if sol(L2) - sol(L1) is larger than some quantity solve the ode for L = L3 > L2 ... and proceed until convergence).
See below.

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

ode:=diff(y(x),x$2)-y(x)=0;

Ls := [$1..6];

for L in Ls do
  IC      := y(0)=1, y(L) = 0;
  sol_||L :=dsolve([ode,IC], numeric);
end do:

col := table(Ls =~ [red, green, blue, black, gold, magenta]):

plots:-display(
  seq(
    plots:-odeplot(sol_||L, [x, y(x)], x=0..L, color=col[L], legend=typeset('L' = L), axis[2]=[mode=log])
    , L in Ls
  )
  , title="Solutions (log scale) for different values of L"
)

diff(diff(y(x), x), x)-y(x) = 0

 

[1, 2, 3, 4, 5, 6]

 

 

S := proc(L, s)
  if s::numeric then
    abs(eval(y(x), (sol_||L)(s)) - evalf(exp(-s)))
  else
    'procname(_passed)'
  end if:
end proc:

plots:-display(
  seq(
    plot(S(L, s), s=0..L, color=col[L], legend=typeset('L' = L), axis[2]=[mode=log])
    , L in Ls
  )
  , title="absolute value of the error compared with the exact solution"
)

 
 

 

Download Numerical_solution.mw

I'm always surprised to see someone writting such a boundary condition, more generally someone considering +∞ as an ordinary, which is not.
A correct writting should be limit(y(x), x=+∞) = y0.

Here is a little discussion 

restart

kernelopts(version)

`Maple 2015.2, APPLE UNIVERSAL OSX, Dec 20 2015, Build ID 1097895`

(1)

ode:=2*x^(1/2)*diff(y(x),x)-y(x) = -sin(x^(1/2))-cos(x^(1/2));

2*x^(1/2)*(diff(y(x), x))-y(x) = -sin(x^(1/2))-cos(x^(1/2))

(2)

# +oo is not a number but a name (see infinity help page), so I replace +infinity by L
# which I will make further to tend to +infinity

ic:=y(L) = y__0;
sol:=dsolve([ode,ic]);
odetest(sol,[ode,ic]);  # regardless the value of L

y(L) = y__0

 

y(x) = cos(x^(1/2))-exp(x^(1/2))*(-y__0+cos(L^(1/2)))/exp(L^(1/2))

 

[0, 0]

(3)

# Now replace x by t*L, make t tend to 1 and make the result tend to +oo

tsol  := eval(rhs(sol), x=t*L):
limit(tsol, t = 1, left);
limit(%, L=+infinity);

y__0

 

y__0

(4)

# Still with x=t*L, make L tend to +oo and make the result tend to +1

limit(tsol, L=+infinity);
limit(%, t = 1, left);    # a continuum of possible limits
 

limit(cos((t*L)^(1/2))-exp((t*L)^(1/2))*(-y__0+cos(L^(1/2)))/exp(L^(1/2)), L = infinity)

 

-2+y__0 .. 2+y__0

(5)

# Thus two different ways to compute the same limit and two different results.
# Consider 'sol' as a function F of x and L.
# What we saw is that the limit of F(x, L) when both x and L tend to +oo, with x < L,
# depends on the path taken to [+oo, +oo], which means F(x, L) doesn't have a limit.
#
# To get an idea of where the previous result came from, see how Maple returns the limit
# of a function whose limit is undefined

limit(cos(x), x=+infinity);

-1 .. 1

(6)

# By the way here is what Maple 2015 returns on your problem

ic:=y(+infinity) = y__0;
sol:=dsolve([ode,ic]):
odetest(sol,[ode,ic]);

y(infinity) = y__0

 

[0, y__0+undefined]

(7)
 

 

Download limit.mw

Note that your implicit question "Shouldn't Maple return no solution?" already holds for  the limit(cos(x), x=+∞) example. In my opinion the -1..1 result is a mathematical nonsense.

What I believe, and even if Maple does accept this syntax, is that you should never write something like y(+infinity) = y0  but adopt a tighter mathematical approach.

@acer  @Elisha

Based on acer's answer here are different ways to improve, IMO, the aesthetic of the figure

(
Remark: ColumnGraph places the first bar at abscissa "0", not "1", meaning acer's picture presents a shift between the bars and their associated tickmarks.
This "0" is not consistent with dataplot or matrixplot conventions (first bar or cell at location "1"). I believe future Maple version should uniformize this.
)

 

different_aesthetics.mw

For instance
       

I'm sure that other improvements can be done using recent Maple's versions (I use Maple 2015)


@Elisha: By the way, transforming a ColumnGraph into a BarGraph, and vice-versa can be done easily using plottools:-transform:

T := textplot(
  [
    [.4, 48, "45%"],   # The ".4" is the half of the default bar width (or height)
    [1.4, 41, "38%"], 
    [2.4, 54, "51%"], 
    [3.4, 70, "67%"], 
    [4.4, 77, "74%"], 
    [5.4, 94, "91%"]
  ], 
  font = ["TIMES", "BOLD", 12]
):

display(
  [
    plottools:-transform((x, y) -> [y, x])(P), 
    T
  ],
  title = "Figure 20: Comparative Effectiveness of Optimal Control Strategies", 
  labels = ["Control Strategies", "Reduction in Coinfection Burden (%)"], 
  labelfont = ["TIMES", "BOLD", 14], 
  titlefont = ["TIMES", "BOLD", 16], 
  axes = boxed, 
  gridlines = true, 
  view = [0..6, 0 .. 100], 
  size = [1000, 650]
)

which produces this

@bashar27

answer updated

@Alfred_F 

Here is a variant where a polygonal domain with non self-intersecting boundary is randomly drawn.
Random_Pick.mw

@Alfred_F 

First of all I missed a colon after restart.
If this does not fix the problem this is a version issue (I use Maple 2015).

If you are interested in Pick's theorem here is a (in my opinion) a very interesting document from the APMEP (Association des Professeurs de Mathématiques de l'Enseignement Public)  APMEP 516

@nm 

which avoids constructing  H_lines and V_lines:

restart:
max_X:=8: 
max_Y:=5:

grid_points:= [seq(seq([n,m],n=0..max_X),m=0..max_Y)]:
the_points := plots:-pointplot(grid_points, symbol=solidcircle, symbolsize=12, color=red):
plots:-display(
  the_points
  , title="my grid"
  , axis[1] = [gridlines = [max_X, color = blue]]
  , axis[2] = [gridlines = [max_Y, color = blue]]
  , scaling=constrained
);

@Ali Guzel  

restart;

with(plots):

T:=60:
Digits:=10:
epsilon:=20:

F:=(x,y)->epsilon*(1-x^2)*y-x:

sys:=diff(X(t),t)=Y(t),diff(Y(t),t)=F(X(t),Y(t)):

vars:={X(t),Y(t)}: ic:=X(0)=0.2,Y(0)=0.1:

pow := [$-4..-2]:

for p in pow do
  sol[p] := dsolve(
              [sys,ic]
              , numeric
              , stepsize=10^p
              , method=classical[foreuler]
              , output=listprocedure
              , maxfun=0
            ):
end do:

col := [red, green, blue]:

display(
  seq(
     odeplot(
       sol[p]
       , [t,X(t)]
       , 0..T
       , axes=normal
       , style=line
       , numpoints=n
       , labels=["t","x"]
       , color=col[abs(p)-1]
       , legend=typeset(h = cat(`#msup(mo("10"),mo("`, p, `"))`))
     )
     , p in pow
  )
)

 

default_sol := dsolve(
                 [sys,ic]
                 , numeric
                 , output=listprocedure
                 , maxfun=0
               ):

display(
  odeplot(
     sol[pow[1]]
     , [t,X(t)]
     , 0..T
     , axes=normal
     , style=line
     , numpoints=n
     , labels=["t","x"]
     , color=col[abs(pow[1])-1]
     , legend=typeset(h = cat(`#msup(mo("10"),mo("`, pow[1], `"))`))
   )
   ,
   odeplot(
     default_sol
     , [t,X(t)]
     , 0..T
     , axes=normal
     , style=line
     , labels=["t","x"]
     , color=gold, thickness=7, transparency=0.7
     , legend="default solver"
   )
  )

 
 

 

Download Are_you_aware_of_convergence_issues.mw

@Alfred_F 

Years ago I bought the first French translation (2015) of this famous book "What is mathematics? An elemenrtary approach to ideas and methods" by Richard Courant and Herbert Robbins (1941).

Its 1st chapter concerns the natural numbers and on line 9 of the the introductory paragraph (page 7) it is written nombres naturels1✻ ... (The footnote  says that the notes tagged with an asterisk are from the book's translator and are thus absent in the vernacular version of the book).
This note (1) can be read page 27 and I reproduce it here (I hope I'm not infringing copyright in doing so).
Note_1.pdf
There is also a discussion about positive, strictly positive and non negative numbers.
 

@Mac Dude 

About the density plot: you could of use the orientation option in the Histogram3D plot command to get this density plot.
But you would have to adjust the lightning model and remove cell boundaries to get a pretty display.
I believed it's better to have a specific procedure dedicated to density plotting.

@emendes 

If it's not too much to ask, could you explain to me what your code is meant to calculate?
(This is more of a curiosity than something that could help me to answer your question.)


 

restart

ic := 1, 2, 3:
rec := rsolve({f(n) = f(n-1) + f(n-2) + f(n-3), seq(f(i)=ic[i], i=1..3)}, {f});

{f(n) = Sum((_R+1)*(1/_R)^n/(3*_R^2+2*_R+1), _R = RootOf(_Z^3+_Z^2+_Z-1))}

(1)

F := unapply(value(rhs(rec[1])), n)

proc (n) options operator, arrow; sum((_R+1)*(1/_R)^n/(3*_R^2+2*_R+1), _R = RootOf(_Z^3+_Z^2+_Z-1)) end proc

(2)

print~([ seq([i, ic[i]], i=1..3), seq([n, F(n)], n=4..10) ]):

[1, 1]

 

[2, 2]

 

[3, 3]

 

[4, 6]

 

[5, 11]

 

[6, 20]

 

[7, 37]

 

[8, 68]

 

[9, 125]

 

[10, 230]

(3)

ic := 'ic';  # unassign the name ic

ic

(4)

tri := proc(ic)
  local rec, F:
  rec := rsolve({f(n) = f(n-1) + f(n-2) + f(n-3), seq(f(i)=ic[i], i=1..3)}, {f});
  F   := unapply(value(rhs(rec[1])), n);
end proc:

ic := [u, v, w];
tri(ic)(n)

[u, v, w]

 

sum(-(-2*_R^2*v+_R^2*w-2*_R*u-_R*v+_R*w+u+v-w)*(1/_R)^n/((3*_R^2+2*_R+1)*_R), _R = RootOf(_Z^3+_Z^2+_Z-1))

(5)

ic := [1, 2, 3];
tri(ic)(n);
map(n -> tri(ic)(n), [$1..10])

[1, 2, 3]

 

sum((_R+1)*(1/_R)^n/(3*_R^2+2*_R+1), _R = RootOf(_Z^3+_Z^2+_Z-1))

 

[1, 2, 3, 6, 11, 20, 37, 68, 125, 230]

(6)

ic := [1$3];
tri(ic)(n);
map(n -> tri(ic)(n), [$1..10])

[1, 1, 1]

 

sum(-(-_R^2-2*_R+1)*(1/_R)^n/((3*_R^2+2*_R+1)*_R), _R = RootOf(_Z^3+_Z^2+_Z-1))

 

[1, 1, 1, 3, 5, 9, 17, 31, 57, 105]

(7)

 


 

Download tripatouille.mw

It would be much simpler to help you if you upload your worksheet (use the big green up-arrow in the menubar, watchout: the file name must not contain special characters).

An advice: verify if the commands you use are all threadsafe (for instance zip is not)

@dharr 

You're right, I have been careless in the conclusion I drew from Descarte's rule of signs.

I got back to what I did to see that in fact I don't use it explicitely (I don't even know why I mentioned it: probably I thought by the time it would help me, but it didn't).
Anyway, I found my error.
At some point I solve a 3r degree equation and Maple, as usual, gives me the roots r1, r2, r3 where r2 and r3 appears in conjugate complex forms, let's say r2 = a+Ib and r3 = a-Ib.
Given these 3 rootsshould be real be real I conclude that b = 0 and thus r2 = r3, so there is a root of multiplicity 2 (for instance y = z).

Replacing z by y in the three symetric functions of (x, y, z) led me to w = 0 and x = 0.
I vaw to quick: in the expression r2 = a+Ib nothing tels me that 'a', or 'b', or both are not imaginary themselves, meaning the shotcut "b must be null" is an unforgiveable error.

I should have taken the time to analyze my work instead of assuming I was on the right track.
As they say, forewarned is forearmed; let’s hope this serves me well as a lesson.

See you next time!

1 2 3 4 5 6 7 Last Page 1 of 42