Doug Meade

 

Doug

---------------------------------------------------------------------
Douglas B. Meade <><
Math, USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Phone: (803) 777-6183 URL: http://www.math.sc.edu

MaplePrimes Activity


These are replies submitted by Doug Meade

Yes, that is one of my ancient sites. The date on the file indicates it is almost 9 years old. The specific MIME types are no longer valid. Here is what we have on our HTTP servers at USC:
application/maple          mw           # Maple (standard) worksheets
application/maple          mws          # Maple (classic) worksheets
application/maplet         maplet       # Maplet worksheets
The note in our mime.types file indicates these were obtained from Maplesoft (as far back as 1996). One fact I learned just this week is that you can have special MIME types for directories of files on the WWW. All you have to do is create a .htaccess file in the directory and include appropriate AddType lines to this file. For example, if you want to use an HTML meta tag to redirect a Maple worksheet to a different location, replace the file corresponding to the old URL with an HTML file with the exact same name (including extension) as existed previously. If this is a classic worksheet, then you would need to put the following in .htaccess in that directory: AddType text/html mws Now, when the server receives a request for a .mws file in this directory the server sends back the current file (with HTML code). Because this content arrives with a MIME type of text/html (not application/maple) it is displayed and processed as an HTML file and the redirection can be completed. If you want to see this in action, visit the URL http://www.math.sc.edu/calclab/MapletsForCalculus/. This should load a "redirection" HTML page and then take you to a new page. (Well, you won't be able to see the new page unless you happen to have a USC login.) If you are interested in the Maplets For Calculus, we do have a page with open access. That URL is http://www.math.sc.edu/calclab/MapletsForCalculus.html. Why do I mention all of this here? Well, maybe William can use the local .htaccess and AddType information to implement a temporary fix to the problem on MaplePrimes. In this case, the .htaccess file will probably look like:
AddType application/maple  mw mws
AddType applicaiton/maplet maplet
I'm sure we'll hear about it if this does not work or if there is a better way to achieve the same result. Doug
David, The basic approach illustrated in the other post works when the basic ODE has an explicit general solution. This is not going to work in the general case when no explicit solution is available in closed form. In such cases it will be necessary to rely solely upon numerical solutions. Here is one way in which this can be accomplished. (Note that this was done in Maple 9.5, but I expect it will work fine in Maple 10 as well.) The basic idea is to create a function (G) that accepts a value for ya, uses this input to compute a numerical solution to y'=f(y,ya,x), y(xa)=ya and returns the value of the constraint function g(y(0),y(xa)). This function can then be evaluated, plotted, and otherwise manipulated just as any other function in Maple. In particular, that means the equation can be solved. Now, an explicit solution is not to be expected, but fsolve can be used to obtain an approximate numerical solution. The only complication to using fsolve is that the function should return unevaluated if the argument is not numeric. (The most complicated part of this implementation is the extraction of the values from the numerical solution. I might be simpler to let dsolve,numeric return a procedure - but this should be marginally slower.) Once the appropriate value for ya is found, the final solution to the IVP can be found and plotted.

> restart;

> with( plots ):

Warning, the name changecoords has been redefined

>

> xa := 1;

> f := (y,ya,x) -> y*ya*x;

> g := (y0,ya) -> sin(y0)*ya+y0*exp(ya)-2;

Maple Equation

Maple Equation

Maple Equation

> G := proc(ya)
  global f, g, xa;
  local ivp, S;
  if not ya::numeric then     # return unevaluated for non-numeric arguments
    return 'procname'(args)
  end if;
  ivp := { diff(y(x),x)=f(y(x),ya,x), y(xa)=ya };
  S := dsolve( ivp, y(x), numeric, output=array([0,xa]) );
  return evalf( g( S[2,1][1,2], S[2,1][2,2] ) );

> end proc:

> #debug(G);

> G(1);

Maple Equation

> G(ya);

Maple Equation

> plot( G(ya), ya=0..1 );

Maple Plot

>

> Ya := fsolve( G(ya)=0, ya );

Maple Equation

> IVP := { diff(y(x),x)=f(y(x),Ya,x), y(xa)=Ya }:

> S := dsolve( IVP, y(x), numeric, output=listprocedure ):

> Y := eval( y(x), S ):

>

> g(Y(0),Y(xa));

Maple Equation

>

> odeplot( S, [x,y(x)], x=0..xa );

Maple Plot

>

This post generated using the online HTML conversion tool
Download the original worksheet

David, The basic approach illustrated in the other post works when the basic ODE has an explicit general solution. This is not going to work in the general case when no explicit solution is available in closed form. In such cases it will be necessary to rely solely upon numerical solutions. Here is one way in which this can be accomplished. (Note that this was done in Maple 9.5, but I expect it will work fine in Maple 10 as well.) The basic idea is to create a function (G) that accepts a value for ya, uses this input to compute a numerical solution to y'=f(y,ya,x), y(xa)=ya and returns the value of the constraint function g(y(0),y(xa)). This function can then be evaluated, plotted, and otherwise manipulated just as any other function in Maple. In particular, that means the equation can be solved. Now, an explicit solution is not to be expected, but fsolve can be used to obtain an approximate numerical solution. The only complication to using fsolve is that the function should return unevaluated if the argument is not numeric. (The most complicated part of this implementation is the extraction of the values from the numerical solution. I might be simpler to let dsolve,numeric return a procedure - but this should be marginally slower.) Once the appropriate value for ya is found, the final solution to the IVP can be found and plotted.

> restart;

> with( plots ):

Warning, the name changecoords has been redefined

>

> xa := 1;

> f := (y,ya,x) -> y*ya*x;

> g := (y0,ya) -> sin(y0)*ya+y0*exp(ya)-2;

Maple Equation

Maple Equation

Maple Equation

> G := proc(ya)
  global f, g, xa;
  local ivp, S;
  if not ya::numeric then     # return unevaluated for non-numeric arguments
    return 'procname'(args)
  end if;
  ivp := { diff(y(x),x)=f(y(x),ya,x), y(xa)=ya };
  S := dsolve( ivp, y(x), numeric, output=array([0,xa]) );
  return evalf( g( S[2,1][1,2], S[2,1][2,2] ) );

> end proc:

> #debug(G);

> G(1);

Maple Equation

> G(ya);

Maple Equation

> plot( G(ya), ya=0..1 );

Maple Plot

>

> Ya := fsolve( G(ya)=0, ya );

Maple Equation

> IVP := { diff(y(x),x)=f(y(x),Ya,x), y(xa)=Ya }:

> S := dsolve( IVP, y(x), numeric, output=listprocedure ):

> Y := eval( y(x), S ):

>

> g(Y(0),Y(xa));

Maple Equation

>

> odeplot( S, [x,y(x)], x=0..xa );

Maple Plot

>

This post generated using the online HTML conversion tool
Download the original worksheet

My mistake. I did not read your question closely. I thought you were looking for the minimum of f(t). Now I see that you are looking for the minimum of f(tau(x)) for x=500..5000. By the chain rule, diff( f(tau(x)), x ) = D(f)(tau(x)) * D(tau)(x) The comments in my original posting were based on the fact that diff( f(t), t ) = (c-s(t))*m/l+x = 0 when s(t) = c+x*l/m My idea was that you would solve this last equation in the same manner that you solved s(t)=c in your second question (to define tau). Now that I understand your problem, what can we say? Looking at the differential equation for f (above) evaluated at t=tau, we see that D(f)(tau(x)) = (c-s(tau(x))*m/l+x = (c-c)*m/l+x = x In particular, this derivative is never zero (when x=500..5000). So, assuming tau is smooth, the minimum occurs either at an endpoint (x=500 or x=5000) or at a point where diff(tau(x),x)=0. You will need to work on this to see how this could be expressed in terms of the original functions in the problem. I'll conclude by showing a function that defines your function tau: > tau := x -> fsolve( S(t)=c, t ); You could then plot tau(x) using > plot( tau, 500..5000 ); Note that I did not write tau(x) or x=500..5000 in the plot command. With my implementation of tau, if I wrote tau(x) Maple would execute the proc and would not be able to solve the equation S(t)=c for a general (unspecified) value of x. There are plenty of other ways to implement tau that avoid this problem - but these require a little more programming that I do not have the time to explain right now. I hope this gives a little better insight into your question. Doug ----------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Dear Work Hard, First, your system involves quite a few unspecified parameters. You need to give numeric values to all parameters before Maple can return a numeric solution to the initial value problem. (These assignments can come before or after the dsolve command, but must be before the output from dsolve is evaluated.) For lack of any better values, I used: > A,rho,theta,r,x,v,c,l,m := $1..9; I think I can answer your questions: > (1) why cannot I evaluate the value of f(t) at t=5 by using f(5)? Because Maple does not have a definition for the function f. What you want to do is the following. First, use output=listprocedure to force Maple to return a list of procedures. (The default is procedurelist - a single procedure that returns a list of information about the solution at the given value of the independent variable.) > solution := dsolve([diff(s(t),t) = A - A * rho^((1 - r^(theta * t)) * x) - v, > diff(f(t),t)=((c - s(t)) / l) * m + x, > diff(h(t),t) = x, > f(0) = 0, s(0) = 0, h(0) = 0], > numeric, > output=listprocedure); To extract the procedures for f (and h and s) use: > F := eval( f(t), solution ); > H := eval( h(t), solution ); > S := eval( s(t), solution ); Now, you can approximate the value for f(5) with: > F(5); > (2) For each fixed x, there are curves s(t), f(t), and h(t). Given > s(tau) = c, I want to find f(tau)=?. How can I do that? Do I need to > find tau first, then find f(tau)? How? With the answer to (1) this is now easy. > tau := fsolve( S(t)=c, t ); > F(tau); > (3) My ultimate goal is for x in a range, say x=500..5000, find the > minimum of f(tau) (since for each x there is a f(tau)). How can I do > that? This requires a little more thought - but not much. I'll not answer this directly for you, but will leave you with these hints. How would you find the minimum if you knew an explicit formula for f? First, you might look at a graph to visually understand the problem. Next you would probably want to find the critical points of f. Look at your differential equations and see how you can obtain an explicit formula for f'(t). Use the ideas in (2) to solve f'(t)=0 (on an appropriate t interval, if necessary). This should get you pretty close to a reasonable answer. I hope this helps, Doug ------------------------------------------------------------------ Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Another way to approach this problem is to use the output= option with an array (or Array) of values for the independent variable. You mention wanting the values of the solution at t=0..100. This can be obtained using:
dsol1 := dsolve({deq1,ic1}, numeric, output=array([$0..100]):
The result will be a 2x1 matrix in which dsol1[1,1] is the column headers, [t, y(t)], and dsol1[2,1] is the 101x3 matrix of values of the solution. Each row of dsol1[2,1] contains three elements: values of t, y(t), and y'(t). That is dsol1[2,1][1,1] willbe 0 (for t=0), dsol1[2,1][1,2] will be the value of y(0), and dsol1[2,1][1,3] will be the value of y'(0). Now, you can use some of the other suggestions to write the numerical solution to a file. I hope this is of some use to you, Doug --------------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Another way to approach this problem is to use the output= option with an array (or Array) of values for the independent variable. You mention wanting the values of the solution at t=0..100. This can be obtained using:
dsol1 := dsolve({deq1,ic1}, numeric, output=array([$0..100]):
The result will be a 2x1 matrix in which dsol1[1,1] is the column headers, [t, y(t)], and dsol1[2,1] is the 101x3 matrix of values of the solution. Each row of dsol1[2,1] contains three elements: values of t, y(t), and y'(t). That is dsol1[2,1][1,1] willbe 0 (for t=0), dsol1[2,1][1,2] will be the value of y(0), and dsol1[2,1][1,3] will be the value of y'(0). Now, you can use some of the other suggestions to write the numerical solution to a file. I hope this is of some use to you, Doug --------------------------------------------------------------------- Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
Just a few weeks ago I discovered the TypeTools package. This package appears to provide enhanced flexibility for users to create their own types within Maple. Here is the code for a procedure to create a new type that tests for a number being in an interval - open or closed, infinite or finite.
> MakeIntervalType := proc(N::symbol, R::range)
>   local Hi, HiEnd, Lo, LoEnd, t, Test, Test1, Test2;
>   description `Creates a user-defined Maple type for an interval - finite or infinite, open or closed`;
>   options `Douglas B. Meade, Univ. of South Carolina, August 2005`;
>   Lo,Hi:=op(R);
>   LoEnd,HiEnd:='closed','closed';
>   if nargs>2 then LoEnd:=args[3] end if;
>   if nargs>3 then HiEnd:=args[4] end if;
>   Test1 := `if`(LoEnd='closed',t>=Lo,t>Lo);
>   Test2 := `if`(HiEnd='closed',t<=Hi,t<Hi);
>   Test := unapply( evalb(Test1 and Test2), t );
>   TypeTools[AddType](N, Test)
> end proc:
Then, to define three different types:
> MakeIntervalType(`(2,8)`,2..8,'open','open');
> MakeIntervalType(`[2,8)`,2..8,'closed','open');
> MakeIntervalType(`(-inf,5]`,-infinity..5,'open');
You can check these definitions as follows:
> TypeTools[GetTypes]();
> TypeTools[GetType](`(2,8)`): showstat(%);
> TypeTools[GetType](`[2,8)`): showstat(%);
> TypeTools[GetType](`(-inf,5]`): showstat(%);
And, finally, to test their functionality:
> for t from -10 to 10 do
>   print( t, type(t,`(2,8)`), type(t,`[2,8)`), type(t,`(-inf,5]`) );
> end do;
I note that MakeIntervalType could be enhanced to incorporate other options, or to automatically generate the type name. I'll leave these enhancements as exercises for the interested reader (but please share your results!) (I stumbled upon TypeTools after I was having trouble using compound types to test contents of a maplet element in a Get command. In that case, TypeTools has worked great.) Prof. Douglas B. Meade Phone: (803) 777-6183 Department of Mathematics URL: http://www.math.sc.edu/~meade/ USC, Columbia, SC 29208 E-mail: mailto:meade@math.sc.edu
First 74 75 76 Page 76 of 76