Carl Love

Carl Love

24830 Reputation

25 Badges

10 years, 115 days
Natick, Massachusetts, United States
My name was formerly Carl Devore.

MaplePrimes Activity

These are answers submitted by Carl Love

See the command ?LinearAlgebra:-GenerateMatrix. After that, you'll want to use ?LinearAlgebra,LinearSolve.

In addition to Adri's suggestions, you could set ?printlevel to a high value. For example

printlevel:= 100;

and then run your code.

Not having your code to work with, my first guess is that examplenr is not properly defined when you read it in another worksheet. I suggest that you construct the filename outside the plotsetup command and look at what it is. Something like this:

filename:= cat(filefolder, examplenr, ".gif");
plotsetup(gif, plotoutput= filename));


When looking at the description of the algorithm in, for example, Wikipedia, you have to think of y (the dependent variable) as a vector; and think of f (the function that evaluates the derivatives y') as a vector-valued function. Try the code below for your procedure Fourth. I kept the structure pretty close to what you had. If I were writing it from scratch, I would explicitly use Vectors.

local f,g,e,n,k1,k2,k3,k4,l2,l1,l3,l4,m1,m2,m3,m4,
   H, X, Y, Z
   f:= (x,y,z)-> a*x-q*x*y-nn*x*z;
   g:= (x,y,z)-> -b*y+k*x*y;
   e:= (x,y,z)-> -c*z+m*x*z;

   H:= h/2;

   for n from 0 to N do
      X:= x[n]; Y:= y[n]; Z:= z[n];
      k1:= f(X,Y,Z); l1:= g(X,Y,Z); m1:= e(X,Y,Z);

      X:= x[n]+H*k1; Y:= y[n]+H*l1; Z:= z[n]+H*m1;
      k2:= f(X,Y,Z); l2:= g(X,Y,Z); m2:= e(X,Y,Z);

      X:= x[n]+H*k2; Y:= y[n]+H*l2; Z:= z[n]+H*m2;
      k3:= f(X,Y,Z); l3:= g(X,Y,Z); m3:= e(X,Y,Z);

      X:= x[n]+h*k3; Y:= y[n]+h*l3; Z:= z[n]+h*m3;
      k4:= f(X,Y,Z); l4:= g(X,Y,Z); m4:= e(X,Y,Z);

      t[n+1]:= t[n] + h;
      x[n+1]:= x[n] + h/6*(k1 + 2*k2 + 2*k3 + k4);
      y[n+1]:= y[n] + h/6*(l1 + 2*l2 + 2*l3 + l4);
      z[n+1]:= z[n] + h/6*(m1 + 2*m2 + 2*m3 + m4)

end proc:

We can use sscanf to break the input string into word-sized chunks, and then pass those words to sprintf so that they can be reassembled with the spaces between words. The syntax of sprintf and printf is identical. The reason that I chose sprintf is that there's no way to programmatically access the output of printf; it produces output for viewing only. 

wordlength:= 5:
# Construct scan format code
WFI:= sprintf("%%%ds", wordlength);


# Construct print format code
WFO:= sprintf("%%0.%ds", wordlength);


# Count words
words:= iquo(length(S), wordlength, 'r'):
if r>0 then  words:= words+1  end if:

sprintf(cat(WFO, cat(" ",WFO) $ words-1), sscanf(S, cat(WFI $ words))[]);

Is this what you were looking for?

The error message about "points with fewer or more than 2 components" clearly provides a clue to the problem. The plot command is two-dimensional (2D), and thus cannot handle the 3D points that you try to give it in the plots axyz and ixyz. For these two plots, change the command plot to pointplot3d.

There's another bug in your code, more subtle and insidious. You use variable n several ways, some of which are in conflict. You use it as a parameter in the first differential equation de1, and as such it also appears in the f:= (x,y,z)-> ... in procedures ESys and Fourth. But you also use it as a local variable (for a loop index) in those procedures. You need to change n to something else, say nn, in four places: the definition of de1, the line where you set the numeric values of the parameters, and the two f:= (x,y,z)-> ... lines.

Consider adding some spaces to your code to make it easier to read.

If you have further questions about this piece of code, please upload a worksheet. This piece is a bit too long for an easy cut-and-paste by the reader (me in this case). After I cut-and-paste, I still have to manually block off the comments and remove the extra >s.

10^5 is not too many points for the plot renderer, but 10^6 is unwieldy. Execute these commands to create the plot:

nn:= nextprime(10^7):
zz:= 1:
N:= 10^5:
Z:= Matrix(N, 2):
for k to N do
   zz:= (zz^2+1) mod nn;
   Z[k,1]:= k:  Z[k,2]:= zz
end do:
plot(Z, style= point);

The plot will appear as a solid block of points. Put the mouse pointer in the plot and right click. Select Manipulator from the context menu, then Scale from the submenu. Select a rectangle on the plot which is horizontally narrow but vertically covers the full extent of the plot. Do this a few more times, until you're satisfied with the narrowness of the horizontal range. Then go to the manipuator submenu again and select Pan. Now the plot is horizontally scrollable.

The problem is that EmpiricalDistribution is a discrete distribution, yet the sample is drawn from a continuous distribution. Something almost analogous to EmpiricalDistribution that works for continuous distributions is KernelDensity. Specifically, KernelDensity returns a probability density function (PDF) deduced from a sample drawn from a continuous distribution. Once you have the PDF, you can create a custom Distribution, and then a RandomVariable from that.



N:= RandomVariable(Normal(10,10)):

P:= Sample(N, 2^9):

A much larger sample would work better, but I'm sticking close to your 500.


F:= KernelDensity(P, method= exact):

plot(F, -20..40);

EmpDist:= Distribution(PDF= F):

X:= RandomVariable(EmpDist):

The DensityPlot is just the same PDF.

DensityPlot(X, range= -20..40);

Compare with the underlying distribution:





Answering your first question, I strongly suspect that in the code before the code that you show that you assigned a plot to the variable a. Doing a restart should fix that.

I believe that the following will work in Maple 12. Let me know if it does not.


a:= 0:  b:= 0.8:  step:= 0.2:

IVP:= {diff(y(t),t)=y(t)-2, y(a)=1}, {y(t)}:

soln1:= dsolve(
   IVP, numeric, method= classical[foreuler],
   stepsize= step,
   output= Array([seq](T, T= a..b, step))

soln1 := Matrix(5, 2, {(1, 1) = 0., (1, 2) = 1., (2, 1) = .2, (2, 2) = .800000000000000, (3, 1) = .4, (3, 2) = .560000000000000, (4, 1) = .6, (4, 2) = .272000000000000, (5, 1) = .8, (5, 2) = -0.735999999999998e-1})

soln2:= dsolve(IVP, numeric):  #high-accuracy numerical solution

soln3:= dsolve(IVP);  #exact solution

y(t) = 2-exp(t)

P1:= plot(
   [soln1 $ 2],
   style= [point, line],
   symbol= cross, symbolsize= 16,
   thickness= 0,
   color= red

P2:= plots:-odeplot(soln2, t= a..b, color= blue):

P3:= plot(eval(y(t), soln3), t= a..b, color= green):





Here (attached) is your code with my (minor) modifications to make it run in hardware floating point. I ran it up to size 512 in less than a minute, and got reasonable accuracy (norm of the residuals 1e-16 seems reasonable to me). Let me know if this fulfills your needs.

In your program, what is the purpose of the following lines?

if type(y, rational) = false then
   y:= evalf(y);
   y:= convert(y, rational)
end if;

In other words, why are you using exact rational arithmetic, even after floating point approximation?

I ran your code for m=50, and the answers have more than 10000 digits in the numerator and denominator. It could be over a million digits for m=300.

Maple's sign is the sign of the leading coefficient of a polynomial, so sign(x) is always 1. Maple's signum is fine if you are dealing with strictly real expressions; for complex values, it may be counterintuitive (I'm not saying that it's wrong---it just may not be what you expect). An alternative is csgn. See ?signum, ?csgn, and ?sign.

Your IVP can be solved symbolically. The solution is x(t) = (1+exp(-t))/2. Thus .5 < x(t) < 1 for all t > 0. Thus both piecewises and the abs make no difference at all and the whole ODE reduces to

diff(x(t), t) = .5 - x(t).

When you make assumptions on a variable, a copy of that variable is made. If you then assign to that same variable, then you'll have the original and the copy being used in different ways, and things get confusing. The takeaway advice for the inexperienced user is Don't assign to assumed variables.

Your last two lines were

  • a:=solve(z1,a);
  • z2;

You should change those to

  • soln:= solve(z1, {a});
  • eval(z2, soln);

When used in this way, eval has most of the desired effect of an assignment statement without the side-effect of changing the global state. You can of course assign the result of eval to another variable.

First 345 346 347 348 349 350 351 Last Page 347 of 363