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

Very nice, but how does Maple TA respond when the student's response is not an initial value problem in correct Maple syntax?

Also, but is checking for equality of solutions sufficient? I would have tried

evalb( simplify( {$sol}-{a} = 0 ) )

Can't you do some more testing of the student's response to provide more instructive feedback that "right" or "wrong"?

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

Your problem appears to be unclear - and incomplete.

Is this your differential equation"

 

To have a hope of uniquely defining the function f(y) you need an initial condition.

It appears as though q is a parameter in this problem, and that A is defined as q*int( f(y), y=-infinity..infinity )

If you have an explicit formula for f and if Maple is able to evaluate the improper integral, then it won't be a problem to create the plot you desire.

But, I'm doubtful that an explicit solution for your problem is possible. (You also need to address the fact that vn and a have not been specified.)

Maple does have some ability to work with parameters in an IVP (see ?dsolve,numeric,interactive). This can work nicely with Maple's numeric solver and the odeplot command (in the plots package).

I think your best hope for creating this plot might be to find a differential equation satisfied by A, as a function of q. If you can do this, then I'm pretty confident you will be able to use odeplot to plot q vs. A.

Here is a very simple problem showing the basic approach I am suggesting:

with( plots ):
IVP := dsolve({diff(y(x), x) = a*y(x), y(x0) = y0}, y(x), numeric, parameters = [x0, y0, a]):
IVP( parameters );
[x0 = undefined, y0 = undefined, a = undefined]
IVP( parameters=[0,1,-1] );    # specific values for the parameters
[x0 = 0., y0 = 1., a = -1.]
IVP(1);
[x = 1., y(x) = HFloat(0.3678793563072187)]
odeplot( IVP, [x,y(x)], x=0..1 );
IVP( parameters=[0,1,-1] );   # another set of parameters
odeplot( IVP, [x,y(x)], x=0..1 );

 I hope some of this is usefull.

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

This problem appears to have very large (10^(17)) and very small (10^(-25)) numbers - all entered as floating-point numbers. The short interval is also an issue. Yet, an analytic solution is desired. 

There does appear to be some structure to the system of DEs. I wonder if it's not possible to work with a general system with the right structure, then insert numerical values for the parameters at the end.

Do you really want a "series" solution, or just the first few terms of a series solution? If it's the latter, maybe you can substitute the desired form of your solution into the system and solve the resulting (linear) system of equations for the unknown coefficients?

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

@vv The curve I found is the relationship between x and y when a=0. (Recall that I used F(0)=10^(-6).) For other values of a, I expect a different relation - I did nothing for these cases. I think there's a lot of information that still remains to be provided to us by the original poster.

@Kitonum Your answer can be explained on geometric grounds, as well. Another description of the set is all points that are closer to 1+I than to -1-I. The boundary of this set is clearly all points on the line where Im(z)=-Re(z).

This result can also be seen in your approach with one additional command:

evalc(eval(abs(z+1+I)>=abs(z-1-I), z=x+I*y));
map(t->t^2, %);

solve( (21), {x,y} );
{y = y, -y <= x}

The requirement that -y<=x is equvalent to y>=-x, with is just Im(y) >= -Re(y), as I claim above.

Doug

@Doug Meade Following up on my previous post:

evalf( [a,b,c] );
[-0.002398966816 - 0.9307744369 I, 0.6665620885 + 1.239067383 I, 1.542028770 - 0.5356027340 I]
evalf( a^2+b^2+c^2 );
0.133664179 + 0.004466829 I

evalf( abs(a)^2+abs(b)^2+abs(c)^2 );
5.510662818

I don't think floating-point errors are involved with this problem.

 

@peacess I don't see a photo in your latest post. Can you upload the actual worksheet?

What do you expect Maple to provide as a response to an input of "X^2+5-2"? It should give "X^2+3" (as below).

X^2+5-2
2
X + 3

Acer

Thank you for the detailed and informative response

I know it seems like a pain to include the ; or : terminator on every command but they do serve a purpose I guess I'm just too much of a programmer and accept this as part of the language (Commas and periods do serve a useful purpose in written English as well)

In this case, the original code came from a MaplePrimes posting I wondered if maybe that was the cause but could not imagine how or why Now I know I thought I was configured for automatic updates but I'll install that update immediately

Thanks again

Doug

PS It's been hard to omit the commas and periods in this response Particularly because I wanted to end with an exclamation point

---------------------------------------------------------------------
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

test2.mw 

Here's the worksheet that I used to generate my original post. It does not have the output showing the error.

I don't see this problem on my setup at home. It is using Maple 2015.2. I'll try it again tomorrow after a reboot.

What is really odd is that when I first discovered this, that worksheet had a proc definition that used local and Maple did not complain. But, it has complained about every use of local afterwards. If I re-execute the proc that accepts the local, Maple does not complain. It's almost as though something in that proc causes this problem. But, restart does not fix the problem - as my example showed. The original worksheet - with output intact is attached next.

test.mw

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

@Annonymouse I'll guess that your main question is how to find the values of the components of your solution at t=500. It would be nice if it was as simple as

sol(500);
[t(500) = 500., x(t)(500) = HFloat(6.515376130777917e217), y(t)(500) = HFloat(1.1415231204727265e109)]

While this is not quite what we might have hoped for, it does provide enough information get what we need.

eval( [x(t)(500),y(t)(500)], sol(500) );
[HFloat(6.515376130777917e217), HFloat(1.1415231204727265e109)]

If this is something that you're going to do often, and possibly with different times and solutions, then you might want a command such as the following:

evalPT := (s,p) -> eval([x(t)(p),y(t)(p)],s(p)):

Now, to plot a point on your solution plot you could do something like this:

sol := dsolve({diff(x(t), t) = piecewise(t <= 500, x(t)+y(t), x(t)-y(t)), diff(y(t), t) = x(t)/y(t), x(0) = 1, y(0) = 1}, {x(t), y(t)}, numeric, output = listprocedure):
P1 := plots:-odeplot(sol, [x(t), y(t)], t = 450 .. 550);
P2 := plots:-pointplot( evalPT(sol,500), symbol=solidcircle, symbolsize=20, color=blue );
plots:-display( P1, P2 );

If you paid attention to the coordinates of the point at t=500 you will see that this point might as well be at the origin on a graph whose axes on the scale of 10^239 (x) and 10^119 (y).

Have you thought about plotting your solution (and point) on a log-log scale?

P3 := plots:-odeplot(sol, [log(x(t)), log(y(t))], t = 0 .. 550);
P4 := plots:-pointplot( log~(evalPT(sol,500)), symbol=solidcircle, symbolsize=20, color=blue );
plots:-display( P3, P4 );

I hope this is helpful.

Doug

@Firehawk120 

Could you tell us what you are computing in your Fortran program? It's very possible that we can do the entire problem in Maple. That would seem to make the most sense, at least until you tell us otherwise.

Doug

@want to be a permanent vegan 

Did you try to use the ideas I wrote about in my original response? If you want to get values of the solution on the interval with steps of size 0.05, you would need to define A := [$0..130*20]/20;. Be warned that list has 2601 entries. You are not going to be able to view this matrix in any meaningful way.

The trouble with I2 is because your eta is a function that returns a (piecewise-defined) function. You probably want to be referring to eta(i), not eta[i], but even this is not going to make sense as the limits of an integral - where are you going to evaluate the function eta(i)?

If, as Carl Love has indicated, all you really want to do is create a plot of the solution, then odeplot is what you should be using. You can control the step size, or let Maple do this for you. I suggest looking at the online help to get some ideas about what can be done with odeplot.

Doug

 The equations you list are not linear. What do you expect the result to be?

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

Acer,

You wrote what I thought. I'm not quite sure how one would avoid this situation. The original setting of the environment variable could be saved, but is there a (good) way to ensure that the correct value is used at the correct times?

Another question is how many other common names are the names of modules that we don't know about? It would probably be enough to adopt a naming convention that used a capital letter as the LAST letter of the name or the second letter of the name. Not the first letter of the name - that's too likely to be used by a programmer.

An alert that alerts users to possibly unexpected consequences to what should be protected names should also be flagged.

Doug

It sounds as though you are asking for the creation of a discont=true option for plot3d.

This would be a nice addition. I can see that it is a difficult feature to implement, but that should not prevent us (and Maplesoft) from working to build the tools needed to make progress in this area.

Until then, plotting points is a classical solution to this problem. And, once you know where the discontinuities are, then you can try to assemble the full picture in pieces. The real question is how much of this can realistically be automated?

Just some thoughts - and a challenge / request.

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
2 3 4 5 6 7 8 Last Page 4 of 76