Robert Israel

6522 Reputation

21 Badges

18 years, 185 days
University of British Columbia
Associate Professor Emeritus
North York, Ontario, Canada

MaplePrimes Activity


These are replies submitted by Robert Israel

One way that may work is to use pade from the numapprox package, which will give you a rational function of specified degrees in numerator and denominator that matches the given terms of the Taylor series.  This will often give a good approximation of the solution over a much larger interval, especially if the singularities in the complex plane are poles.  In this case:

> F:= (-1)-1+x+(1/2)*(-1+x)^2+(-1+x)^3+(2/3)*(-1+x)^4+(37/60)*(-1+x)^5+(127/360)*(-1+x)^6+(31/140)*(-1+x)^7
 +(185/2016)*(-1+x)^8+(53/1620)*(-1+x)^9-(391/151200)*(-1+x)^10-(1097/90720)*(-1+x)^11
 -(160903/11975040)*(-1+x)^12-(969289/97297200)*(-1+x)^13-(4366883/681080400)*(-1+x)^14
 -(5242229/1513512000)*(-1+x)^15;
  P87:= numapprox[pade](F,x=1,[8,7]);
  plot(P87, x=0 .. 4);

One way that may work is to use pade from the numapprox package, which will give you a rational function of specified degrees in numerator and denominator that matches the given terms of the Taylor series.  This will often give a good approximation of the solution over a much larger interval, especially if the singularities in the complex plane are poles.  In this case:

> F:= (-1)-1+x+(1/2)*(-1+x)^2+(-1+x)^3+(2/3)*(-1+x)^4+(37/60)*(-1+x)^5+(127/360)*(-1+x)^6+(31/140)*(-1+x)^7
 +(185/2016)*(-1+x)^8+(53/1620)*(-1+x)^9-(391/151200)*(-1+x)^10-(1097/90720)*(-1+x)^11
 -(160903/11975040)*(-1+x)^12-(969289/97297200)*(-1+x)^13-(4366883/681080400)*(-1+x)^14
 -(5242229/1513512000)*(-1+x)^15;
  P87:= numapprox[pade](F,x=1,[8,7]);
  plot(P87, x=0 .. 4);

> with(plots):
  animate(plot3d, [max(0, t - sqrt((x-1)^2 + (y-2)^2), t - sqrt((x-3)^2 + (y-3)^2), t - sqrt((x-4)^2 + (y-1)^2)), 
     x = -4 .. 8, y = -4 .. 8, style=patchnogrid, shading=zhue], t = 0 .. 2);
> with(plots):
  animate(plot3d, [max(0, t - sqrt((x-1)^2 + (y-2)^2), t - sqrt((x-3)^2 + (y-3)^2), t - sqrt((x-4)^2 + (y-1)^2)), 
     x = -4 .. 8, y = -4 .. 8, style=patchnogrid, shading=zhue], t = 0 .. 2);

?odeplot works fine for me in Maple 13.02 and Maple 14.  In Maple 12 Classic you need ?plots,odeplot. 

But odeplot is for an output of dsolve(..., numeric), not dsolve(..., series).

?odeplot works fine for me in Maple 13.02 and Maple 14.  In Maple 12 Classic you need ?plots,odeplot. 

But odeplot is for an output of dsolve(..., numeric), not dsolve(..., series).

Or:

min(select(type, [S1,S2], numeric));

This should be a somewhat more robust solution, as it should discard any sort of non-numeric items, not just FAIL

Or:

min(select(type, [S1,S2], numeric));

This should be a somewhat more robust solution, as it should discard any sort of non-numeric items, not just FAIL

If the problem doesn't supply initial and boundary conditions, you might try the symbolic version of pdsolve, which will give you some solutions (not a complete set in any way):

> sys:= {diff(A,t)^2 - N^2*diff(A,y$2) - N^2*diff(A,y)^2 = 0,
   N^3*diff(A,y)^2 + N^2*A*diff(A,y)*diff(N,y) -N*A*diff(A,t$2)- N*A*diff(A,t)^2+A*diff(A,t)*diff(N,t)=0};
  pdsolve(sys);

{a(t,y) = 0, n(t,y) = n(t,y)}, {a(t,y) = _C1, n(t,y) = n(t,y)}, {a(t,y) = ln(_C1*y+_C2), n(t,y) = 0}, {a(t,y) = ln(_C1*y+_C2), 
  n(t,y) = _F1(t)/ln(_C1*y+_C2)}, {a(t,y) = _F1(y), n(t,y) = 0}

You can also try casesplit in the PDEtools package, which will try to "decouple" the equations.

> C:= PDEtools[casesplit](sys);

 The first case is probably what you want: it has n^2 expressed in terms of  derivatives of a (from your first equation), and then a rather complicated pde for a(t,y).

 > Neq:= op([1,1],C[1]);

Neq := n(t,y)^2 = diff(a(t,y),t)^2/(diff(a(t,y),`$`(y,2))+diff(a(t,y),y)^2)

> Aeq:= op([1,2],C[1]);

Aeq := diff(a(t,y),t,`$`(y,2)) = (-2*a(t,y)*diff(a(t,y),t)*diff(a(t,y),`$`(y,2))^2-6*diff(a(t,y),t)*diff(a(t,y),`$`(y,2))*diff(a(t,y),y)^2*a(t,y)-2*diff(a(t,y),t)*diff(a(t,y),y)^4*a(t,y)+2*diff(a(t,y),t)*diff(a(t,y),y)^4+2*diff(a(t,y),t)*diff(a(t,y),`$`(y,2))*diff(a(t,y),y)^2-a(t,y)*diff(a(t,y),t)*diff(a(t,y),y)*diff(a(t,y),`$`(y,3)))/(a(t,y)*diff(a(t,y),y)^2+a(t,y)*diff(a(t,y),`$`(y,2)))

Now you can try solving this pde numerically for a(t,y) with appropriate initial and boundary conditions.  Note that it is first order in t, and third order in y, so one initial and three boundary conditions.

> ibc1:= {a(0,y)=1,a(t,0)=1,D[2](a)(t,0)=1,a(t,1)=t};
  Sol1:= pdsolve(Aeq,ibc1,numeric);

Unfortunately, this equation may be hard to solve.

> Sol1:-plot3d(A,t=0..1,y=0..1,axes=box);

Error, (in pdsolve/numeric/plot3d) unable to compute solution for t>0.:
solution becomes undefined, problem may be ill posed or method may be ill suited to solution

But I was able to find some initial and boundary conditions that worked.

> ibc2:= {a(0,y)=y,a(t,0)=0,D[2](a)(t,0)=1+t,a(t,1)=1+2*t};
  Sol2:= pdsolve(Aeq,ibc2,numeric);
  Sol2:-plot3d(A,t=0..1,y=0..1,axes=box);

If the problem doesn't supply initial and boundary conditions, you might try the symbolic version of pdsolve, which will give you some solutions (not a complete set in any way):

> sys:= {diff(A,t)^2 - N^2*diff(A,y$2) - N^2*diff(A,y)^2 = 0,
   N^3*diff(A,y)^2 + N^2*A*diff(A,y)*diff(N,y) -N*A*diff(A,t$2)- N*A*diff(A,t)^2+A*diff(A,t)*diff(N,t)=0};
  pdsolve(sys);

{a(t,y) = 0, n(t,y) = n(t,y)}, {a(t,y) = _C1, n(t,y) = n(t,y)}, {a(t,y) = ln(_C1*y+_C2), n(t,y) = 0}, {a(t,y) = ln(_C1*y+_C2), 
  n(t,y) = _F1(t)/ln(_C1*y+_C2)}, {a(t,y) = _F1(y), n(t,y) = 0}

You can also try casesplit in the PDEtools package, which will try to "decouple" the equations.

> C:= PDEtools[casesplit](sys);

 The first case is probably what you want: it has n^2 expressed in terms of  derivatives of a (from your first equation), and then a rather complicated pde for a(t,y).

 > Neq:= op([1,1],C[1]);

Neq := n(t,y)^2 = diff(a(t,y),t)^2/(diff(a(t,y),`$`(y,2))+diff(a(t,y),y)^2)

> Aeq:= op([1,2],C[1]);

Aeq := diff(a(t,y),t,`$`(y,2)) = (-2*a(t,y)*diff(a(t,y),t)*diff(a(t,y),`$`(y,2))^2-6*diff(a(t,y),t)*diff(a(t,y),`$`(y,2))*diff(a(t,y),y)^2*a(t,y)-2*diff(a(t,y),t)*diff(a(t,y),y)^4*a(t,y)+2*diff(a(t,y),t)*diff(a(t,y),y)^4+2*diff(a(t,y),t)*diff(a(t,y),`$`(y,2))*diff(a(t,y),y)^2-a(t,y)*diff(a(t,y),t)*diff(a(t,y),y)*diff(a(t,y),`$`(y,3)))/(a(t,y)*diff(a(t,y),y)^2+a(t,y)*diff(a(t,y),`$`(y,2)))

Now you can try solving this pde numerically for a(t,y) with appropriate initial and boundary conditions.  Note that it is first order in t, and third order in y, so one initial and three boundary conditions.

> ibc1:= {a(0,y)=1,a(t,0)=1,D[2](a)(t,0)=1,a(t,1)=t};
  Sol1:= pdsolve(Aeq,ibc1,numeric);

Unfortunately, this equation may be hard to solve.

> Sol1:-plot3d(A,t=0..1,y=0..1,axes=box);

Error, (in pdsolve/numeric/plot3d) unable to compute solution for t>0.:
solution becomes undefined, problem may be ill posed or method may be ill suited to solution

But I was able to find some initial and boundary conditions that worked.

> ibc2:= {a(0,y)=y,a(t,0)=0,D[2](a)(t,0)=1+t,a(t,1)=1+2*t};
  Sol2:= pdsolve(Aeq,ibc2,numeric);
  Sol2:-plot3d(A,t=0..1,y=0..1,axes=box);

If I'm not mistaken, you already have three initial conditions: a(0,y) = 0,  D[1](a)(0,y) = 1, and n(0,y) = 1, so you don't need another.  You also said you had one equation second order in y and another first order in y, so you need 3 boundary conditions too.  In your second posting in this thread you mentioned three: a(t,0)=1, a(t,1)=2,  n(t,0)=1, so that should be enough.

But I'm curious: where does this problem come from?  Is there something you are modeling (in which case the initial and boundary conditions should come from that something)?


If I'm not mistaken, you already have three initial conditions: a(0,y) = 0,  D[1](a)(0,y) = 1, and n(0,y) = 1, so you don't need another.  You also said you had one equation second order in y and another first order in y, so you need 3 boundary conditions too.  In your second posting in this thread you mentioned three: a(t,0)=1, a(t,1)=2,  n(t,0)=1, so that should be enough.

But I'm curious: where does this problem come from?  Is there something you are modeling (in which case the initial and boundary conditions should come from that something)?


Well, the idea is that we want a parametric curve x = X(t), y = Y(t) with [X(j), Y(j)] approximating a[j] for j from 1 to n.  Since a[1] = a[n], it seems reasonable to make X(t) and Y(t) periodic with period n-1.  The natural thing to look at when talking about periodic functions is trigonometric polynomials:
sum(A[j]*cos(2*Pi*i*j/(n-1)) + B[j]*sin(2*Pi*i*j/(n-1)), j=0..k).  I rather arbitrarily took k= 4, and used Fit to get best least-squares fits to both x and y coordinates.

Well, the idea is that we want a parametric curve x = X(t), y = Y(t) with [X(j), Y(j)] approximating a[j] for j from 1 to n.  Since a[1] = a[n], it seems reasonable to make X(t) and Y(t) periodic with period n-1.  The natural thing to look at when talking about periodic functions is trigonometric polynomials:
sum(A[j]*cos(2*Pi*i*j/(n-1)) + B[j]*sin(2*Pi*i*j/(n-1)), j=0..k).  I rather arbitrarily took k= 4, and used Fit to get best least-squares fits to both x and y coordinates.

Hmm, now you have an end condition as well as initial conditions.  pdsolve(..., numeric) is unable to handle those.  There may not be
any way to do this in Maple, unless you want to program your own numerical method for solving this.

First 32 33 34 35 36 37 38 Last Page 34 of 187