/* Maxima is a GNU free software package for algebraic computations.	  */
/* For downloads and information go to:	 	*/
/*        http://maxima-project.org/		*/
/* The following Maxima code was developed by Tyler Smith. 	*/
/* If using the graphical interface wxMaxima, replace the plot commands below with wxplot. */


/* Let's first practise using the following commands:

    factor[f]  - Factors the function f
    diff[f, x]    - The derivative of a function f with respect to x
    integrate[f, x]  - The indefinite integral of f with respect to x
    integrate[f, x, a, b]  - The definite integral of f with respect to x between a and b

For example, let's start by factoring (1 + x - 6 * x^2)/(1 - 2 * x):  */

factor((1 + x - 6 * x^2)/(1 - 2 * x));

/*  In this equation, asterisks are used to indicate multiplication and ^ is used to indicate a power. Maxima should
  answer with 1 + 3x. Use % as a shortcut for the previous result. Similarly,  %on and %in represent the nth output and
  input respectively:  */

/* factor the first output value */
factor(%o1 + 2 * x^2); 

/* expand previous result */
expand(%); 

	   /* Derivatives: */

/* The first derivative of the expression wrt x */
diff(4 * x^3, x); 

/* assign the result to variable tmp */
tmp: diff(exp(4 * x^3), x); 

/* evaluate previous result for x == 1 */
ev(%, x=1); 

/* previous result as a real number */
float(%); 

/* integration wrt x */
integrate(tmp, x); 

/* just display the integration */
'integrate(tmp, x); 

/*definite integral from x == 0 to x == 1*/
integrate(1 + x^2, x, 0, 1);
 

	/* Solutions to practise questions: */

factor( 1 + 3 * x + 2 * x^2);
factor( 2 - 5 * x - 12 * x^2);
diff(1 + 3 * x + 2 * x^2, x);
diff(x^3 * 3 * y^2, x);
diff(x^3 * 3 * y(x)^2, x);
integrate( 1 + 3 * x + 2 * x^2, x);
integrate(x^3 * 3 * y^2, x);
integrate(x^3 * 3 * y(x)^2, x);

integrate( 1 + 3 * x + 2 * x^2, x, 0, 1);
integrate( (1 - x^2)/(1 - x) , x, 0, 1);
integrate(x^3 * 3 * y^2, x, 0, 1);
integrate(4 * x^3, x, 0, 1);

	    /* Finding the right commands */
/* find an exact match */
? term

/* find a fuzzy match */
?? term

/* basic plotting functions */
plot2d(1 + 3 * x + 2 * x^2, [x, 0, 100]);
plot3d(1 + 3 * x + 2 * x * y, [x, 0, 100], [y, 0, 100]);

/* random integer from 1 to 100: */
random_discrete_uniform(100);


/* random float from 0 to 1: */
random_continuous_uniform(0,1);

/* add a third argument to generate a list of random numbers: */
random_discrete_uniform(100, 10);
random_continuous_uniform(0, 1, 10);

/* alternately, random integer from 0 to 99 (only one at a time though): */
random(100);

	/* Lists: */
mylist: [1, 2, 3, 5, 7, 11];

/* make a list from the sequence 1/n for n = 1 thru 10: */
mynewlist: makelist (1/i, i, 1, 10);

/* no explicit listplot function exists: */
xvals: makelist(i, i, 1, length(mynewlist));

plot2d([discrete, xvals, mynewlist], [style, points]);
plot2d([discrete, xvals, mynewlist]);

/* make our own listplot: */
listplot(x) := (xl: length(x), xval: makelist(i, i, 1, xl), plot2d([discrete, xvals, x], [style, points]));

/* 50 random normal variables, mean = 0 , sd = 1 */
Y: random_normal(0, 1, 50);
X: makelist(i, i, 1, length(Y));
rantab: transpose(matrix(X, Y));

/* 50 random integers, normal with mean 0 and sd 1 */
rantab: apply('matrix, create_list([i, random_normal(0,1)], i, 1, 100));
lsquares_estimates(rantab, [x, y], y = a*x + b, [a, b]); 

/* here's the example from mbe: */
/* build a line with random noise, save it to file: */

doline() := block ([x, y, r, display2d:false], set_random_state(make_random_state (654321)),
 	  with_stdout("fit1.dat", for i : 1 thru 10 do (
	  			  x:i,
				  r: -0.2 + random(0.4),
				  y: 1 + x + r,
				  print(x,y))))$

doline()$
datamatrix: read_matrix("fit1.dat");
load(lsquares)$
lsquares_estimates( datamatrix, [x, y], y = a*x + b, [a, b]);


sinplus (x) := sin(x) + x;


/** solving recursion equations: **/
load("solve_rec");
solve_rec(n[t]=a*n[t-1], n[t], n[0]=n0);
solve_rec(n[t]=n[t-1] + r * n[t-1] * ( 1 - n[t-1]/K), n[t], n[0]=n0);

/** solving differential equations: **/
load("dynamics");
eq: 'diff(n(t), t) = r * n(t) * (1 - n(t)/K);
rk(eq, n, 10, [t, 0, 30, 1]);
desolve('diff(f(n), n, 1) = r * n * (1 - n/K), f(n));

/* To quit: */
quit();