/*	PRUEBA EVALUABLE DE PROGRAMACION EN WXMAXIMA, convocatoria de junio de 2016
	Física Computacional I, 1er curso del grado en física
	Dept. de Física Matemática y de Fluidos, UNED, curso 2015-2016

	Este archivo contiene el código en Maxima de las funciones pedidas en la prueba evaluable.
	Para escribirlo se ha empleado el editor de textos "KWrite" en un PC con Linux Fedora 23.

	Para cargar el archivo desde wxMaxima se puede ejecutar: 

		batchload( concat( path, "2016-J-R.mc" ) );

	donde la variable path es una cadena que indica la localización del archivo
	2016-J-R.mc en el disco, p. ej. path podría ser algo parecido a 

		path : "/home/usuario/Fisica-Computacional-1/Maxima/examen/";
*/

/*	EJERCICIO 1
 
	Input:

		intervalo	=	intervalo [a, b] donde está definido el problema.
		N		=	Número de valores de x considerado en la discretización.
		CC		=	Lista de condiciones de contorno [fa, fb].
		f0		=	Función f0(t) empleada para generar la lista de condiciones iniciales.

	Output:

		[xlist, flist, fxxlist, f0list]

		xlist		=	Lista de valores discretizados de x excluyendo los extremos "a" y "b".
		flist		=	Lista da funciones fi = f(xi, t).
		fxxlist		=	Lista de lados derechos de las ecuaciones de evolución dfi/dt = ...
		f0list		=	Lista de valores iniciales fi(0) = f0(xi).
					Aunque el enunciado no pide generar esta lista aquí, dado que hace 
					falta como input de la función del Ejercicio 2 se ha optado por
					generar en esta función el listado de valores fi(0). Por supuesto 
					que también es válido generar estos valores en el Ejercicio 2 o en 
					otra función independiente.

		OBS:
			En el output de las listas flist[i], fxxlist[i] y f0list[i] vamos a considerar
			explusivamente valores del índice i entre i=2 e i=N-1, ya que los valores extremos 
			(i = 1 e i = N) no se necesitan en el ejercicio 2. Análogamente, para la lista xlist[i]
			esta función devuelve solamente los valores de xi correspondientes a las fi anteriores, 
			es decir eliminamos de la lista de valores xi los extremos x[1]=a y x[N]=b.
*/

funcion1( intervalo, N, CC, f0 ) := block( [h, xlist, flist, fxxlist, f0list],

	/*
		Definimos la lista de valores discretizados de la variable x, según se indica en el enunciado
	*/

	h : (intervalo[2] - intervalo[1]) / (N - 1),

	xlist : makelist( intervalo[1] + h * (i - 1), i, 1, N, 1 ),

	/*
		Definimos la lista de funciones fi(t) según se ha indicado en el foro de consultas de la asignatura
	*/

	flist : makelist( concat(f, i), i, 1, N, 1 ),

	/*
		Asignamos los valores de f1(t) y fN(t) a las condiciones de contorno del problema (fa y fb),
		según indica el enunciado
	*/

	flist[1] : CC[1],

	flist[N] : CC[2],

	/*
		Generamos la lista de valores aproximados de d^2f/dx^2|_x=xi (i = 2, ..., N-1), 
		según indica el enunciado
	*/

	fxxlist : makelist( flist[i+1] - 2*flist[i] + flist[i-1], i, 2, N-1, 1 ) / h^2,

	/*
		Al haber asignado previamente flist[1]=fa y flist[N]=fb la anterior lista de valores de d^2f/dx^2|_x=xi
		ya incorpora las condiciones de contorno del problema.
	*/

	/*
		Generamos la lista de valores iniciales fi(0) = f0(xi) (i = 2, ..., N-1), según indica el enunciado

		OBS:

			El enunciado del problema asume que operamos con una ecuación cuyas variables independientes son
			"x" y "t", siendo "x" la variable espacial y "t" la variable temporal.
			Consideramos por tanto como "variables globales" a dichas variables y en la instrucción siguiente
			asumimos que hemos denominado "x" a la variable de tipo "espacio" y suponemos que f0 es una 
			función (o una expresión) dependiente de la variable "x". 
			La instrucción siguiente no funcionará en problemas donde la variable de tipo espacio no se llame 
			"x", en ese caso bastaría con añadir a los argumentos de esta función un argumento adicional con 
			el nombre de dicha variable.
	*/

	f0list : makelist( subst(xlist[i], x, f0), i, 2, N-1, 1 ),

	/*
		Según hemos indicado más arriba, eliminamos de las listas xlist y flist los valores extremos empleando
		para ello la función delete.
		Organizamos el output de esta función como una lista que contiene todas las listas pedidas:

			[xlist, flist, fxxlist, f0list]
	*/

	xlist : delete( last(xlist), delete( first(xlist), xlist ) ),

	flist : delete( last(flist), delete( first(flist), flist ) ),

	return( [xlist, flist, fxxlist, f0list] )

)$

/*	EJERCICIO 2

	Input:

		flist		=	Lista da funciones fi = f(xi, t).
		fxxlist		=	Lista de lados derechos de las ecuaciones de evolución dfi/dt = ...
		f0list		=	Lista de valores iniciales fi(0) = f0(xi).
		varindep	=	t.
		intervalo	=	intervalo [t inicial, t final] de valores de t donde está definido el problema.
		pasork		=	Paso de integración para el algoritmo de RK.
	Output:

		Matriz con los resultados numéricos obtenidos por medio de la función rk.
*/

funcion2( flist, fxxlist, f0list, varindep, intervalo, pasork ) := block( [ aux ],

	/*	Resolvemos el sistema por medio de rk con el paso indicado en pasork y guardamos la solución en la 
		variable local aux.
		Dado que en los apuntes y PECs de otros cursos ya hay multitud de ejemplos de uso de la función rk,
		no consideramos necesario explicar el significado de cada uno de los argumentos de rk.
	*/

		load(dynamics),

		aux : rk(fxxlist, flist, f0list, [varindep, intervalo[1], intervalo[2], pasork]),

	/*	La solución numérica que nos proporciona rk es una lista de listas de valores.
		Es conveniente convertir el output de rk a formato matricial, para ello convertimos
		esta lista en una matriz (ver soluciones tema 2.10 y examen 2011).
	*/

		aux : transpose( apply(matrix, aux) ),

		return(aux)
)$

/*	EJERCICIO 3

	Tal y como se ha aclarado en el foro de la asignatura, en esta función no es necesario suministrar los valores
	mínimo y máximo de la variable independiente t, ya que dichos valores están implícitos en el output de la 
	función del ejercicio anterior.

	Input:

		xlist	=	Lista de valores x_i generamos en el ejercicio 1.
		outf2	=	Output de la función del ejercicio 2.

	Output:

		Gráfica 3D de f(x, t) de acuerdo a la solución aproximada obtenida por medio de rk.
*/


funcion3( xlist, outf2 ) := block( [ tlist ],

	/*
		Definimos: tlist = lista da valores discretos de t empleados por rk
	*/

	tlist : outf2[1],

	/*
	OBS:
		Realizamos una representación 3D de los puntos fi(xi, tj) de acuerdo a las instrucciones aportadas
		en el foro de la asignatura, por medio del paquete "draw", empleando la función "points" y draw3d".

		De todas las opciones ofrecidas en la ayuda de wxmaxima vamos a emplear la función points con la sintaxis

			points(1d_x_array, 1d_y_array, 1d_z_array)

		donde

			1d_x_array = lista de valores de x

			1d_y_array = lista de valores de t

			1d_z_array = lista de valores de f(x, t)

		Para cada uno de los valores de x considerados (xlist[i], con i fijo) vamos a generar el conjunto de 
		puntos a representar (x_i, t_j, f_ij), recorriendo todos los valores de t_j generadosr por rk, haciendo:

			points(

				makelist( xlist[i], j, 1, length(tlist), 1 ),	= valor x_i constante para cada f_i

				tlist,						= lista de valores de tiempo t_j

				outf2[1+i]					= lista de valores de f_i(t_j) = f(x_i, t_j)

			)

		Repetimos esta instrucción por medio de un bucle, incrementando i en una unidad en cada iteración.
		Finalmente empleamos apply( draw3d, points(...) ) para aplicar la función draw3d sobre la colección de
		puntos generada en el bucle, obteniendo así la gráfica pedida.

		La gráfica obtenida de esta forma muestra puntos discretos (x_i, t_j, f(x_i, t_j) = f_i(t_j)), de todas
		formas, dado que el espaciado entre valores de t suele ser mucho menor que el espaciado entre valores de 
		x, el resultado da la impresión de ser una colección de líneas (x_i, t, f_i(t)).

		De acuerdo a como se ha definido la función 1, al eliminarse los extremos (i=1 e i=N) en las listas de 
		valores de x_i y f_i, en esta función vamos a representar solamente los valores de f correspondientes a
		valores de x interiores al intervalo [a, b], es decir, los valores extremos f(a, t) y f(b, t) (dados 
		por las condiciones de contorno del problema fa(t) y fb(t)) no se representan. Si se quieren representar
		dichos valores basta con añadir las líneas (a, t, fa(t)) y (b, t, fb(t)).
	*/

	load(draw),

	apply( draw3d,

	       makelist(

			points(

				makelist( xlist[i], j, 1, length(tlist), 1 ),

				tlist,

				sol[1+i]

			)

		, i, 1, length(xlist) )

	)

)$

/*	EJERCICIO 4

	Input:

		xlist	=	Lista de valores x_i generamos en el ejercicio 1.
		outf2	=	Output de la función del ejercicio 2.
		K	=	Número de valores de t mayores que t=tinicial a representar.

	Output:

		Gráfica 2D de f(x, t_i) vs x K+1 curvas, correspondientes a los K valores de t mayores que el valor inicial junto con la condición inicial (f0(x)), de acuerdo a la solución aproximada obtenida por medio de rk.
*/

funcion4( xlist, outf2, K ) := block( [ tlist, numvalorest, aux_expo, indicesEj4, data ],

	/*
		En esta función también son aplicables los comentarios realizados para la función del ejercicio 3.
		El único input realmente necesario es la lista de valores de x y el output de la función rk.
		Dado el valor K, suministrado como argumento a esta función, representaremos 1+K curvas f(x, tj),
		correspondiendo la primera de ellas a la condición inicial (f(x, 0) = f0(x)) y las K curvas restantes
		a K valores de t entre t inicial y t final.

		Definimos: tlist = lista da valores discretos de t empleados por rk
	*/

	tlist : outf2[1],

	numvalorest : length(tlist),

	/*
		Tal y como se ha explicado en los foros de la asignatura, introducimos como variable local el exponente
		auxiliar (siempre positivo)

			aux_expo

		que emplearemos para definir la lista de valores discretos de t de forma que t^(1/aux_expo) está equiespaciado entre 0 y tfin.
		De esta forma, si aux_expo > 1 (respectivamente < 1) los valores de t empleados para generar esta familia 
		de curvas están distribuidos de manera más densa cerca de valores pequeños (respectivamente grandes) de t.
		Si tomamos aux_expo = 1 obtenemos una lista de 1+K valores de t equiespaciados entre 0 y tfin.

		Aunque esta forma de mostrar los resultados es interesante, dado que el enunciado no especifica nada al
		respecto también consideraremos como válida cualquier otra forma alternativa de tomar los K valores de t 
		a representar (o 1+K si incluimos también la condición inicial). En este sentido la opción más sencilla es
		la aportada en los foros de la asignatura:

			indicesEj4 : makelist( i, i, 1, numvalorest, round(numvalorest/K) )

		en cuyo caso visualizaremos 1+K curvas f(x, tj) vs x para 1+K valores tj equiespaciados entre la condición
		inicial (incluida) y el valor final de t.
	*/

	aux_expo : 1,

	indicesEj4 : makelist(

			round( 1 + numvalorest * (i/numvalorest)^aux_expo ),

		i, 0, numvalorest - 1, round(numvalorest/K) ),

	data : makelist(

			[ discrete, xlist,

				makelist( outf2[ 1+i, indicesEj4[j] ], i, 1, length(xlist), 1 )

			]

		, j, 1, length(indicesEj4), 1 ),

	/*
		Si se quiere visualizar la lista de tj considerados para las graficas de esta función basta con evaluar 

			tlist : makelist( tlist[ indicesEj4[j] ], j, 1, length(indicesEj4), 1 ),

		Aunque el enunciado no lo pedía, en esta función vamos a generar una leyenda en la gráfica con los valores
		de los tiempos tj visualizados:
	*/

	tlist : makelist( tlist[ indicesEj4[j] ], j, 1, length(indicesEj4), 1 ),

	plot2d(data,

		append( [legend], makelist(

				concat( "t = ", string(tlist[i]) ),

			i, 1, length(tlist), 1 ) )

	)

)$

/*
	Como trabajo adicinal definimos la función 4 con extremos funcion4masEXT exactamente igual que la función
	anterior, pero recibiendo como argumentos extra el intervalo [a, b] de definición de x y las condiciones de 
	contorno del problema, lo cual nos permite representar también los valores correspondiendo a los extremos del 
	intervalo.

	El input de esta función es como el de la función 4, con los argumentos adicionales intervalo y CC:

		intervalo	=	intervalo [a, b] donde está definido el problema.
		CC		=	Lista de condiciones de contorno [fa, fb].

	Por claridad en el código eliminamos de esta función los comentarios realizados en la función 4.
*/

funcion4masEXT( xlist, outf2, K, intervalo, CC ) := block( [ tlist, numvalorest, aux_expo, indicesEj4, data ],

	tlist : outf2[1],

	numvalorest : length(tlist),

	aux_expo : 1,

	indicesEj4 : makelist(

			round( 1 + numvalorest * (i/numvalorest)^aux_expo ),

		i, 0, numvalorest - 1, round(numvalorest/K) ),

	tlist : makelist( tlist[ indicesEj4[j] ], j, 1, length(indicesEj4), 1 ),

	data : makelist(

			[ discrete,

				/*  añadimos a xlist los valores inicial y final del intervalo  */

				append( [intervalo[1]], xlist, [intervalo[2]] ),

				/*  añadimos a fi(tj) los valores fa(tj) y fb(tj)  */

				append( 

					[ subst( tlist[j], t, CC[1] ) ],

					makelist( outf2[ 1+i, indicesEj4[j] ], i, 1, length(xlist), 1 ),

					[ subst( tlist[j], t, CC[2] ) ]

				)

			]

		, j, 1, length(indicesEj4), 1 ),

	plot2d(data, append( [legend], makelist( concat( "t = ", string(tlist[i]) ), i, 1, length(tlist), 1 ) ) )

)$

/*
	Para ver cómo funcionan estas funciones se puede ejecutar las instrucciones siguientes en una sesión de wxMaxima:

		intervalo_x : [0, 1];

		intervalo_t : [0, 10];

		CC : [t/(1+t), 0];

		fini : 0;

		N : 25;

		evpaso : 0.001;

		of1 : funcion1( intervalo_x, N, CC, fini )$

		sol : funcion2( of1[2], of1[3], of1[4], t, intervalo_t, evpaso )$

		funcion3( of1[1], sol );

		funcion4( of1[1], sol, 20 );

		funcion4masEXT( of1[1], sol, 20, intervalo_x, CC );

	Este problema describe el calentamiento inestacionario de una barra de longitud L, inicialmente a temperatura 
	T = T1, en la que uno de sus extremos se calienta de acuerdo a la función 

		T = T1 + (T2 - T1) * t / (1 + t)

	mientras que el otro extremo se mantiene a la temperatura inicial T1.

	El problema se resuelve en términos de la coordenada adimensional, x, definida como 

		x = distancia recorrida desde el extremo 1 de la barra / L

	y la temperatura adimensional 

		f = (T - T1) / (T2 - T1)

	A tiempos cortos la temperatura del extremo 1 va aumentando de acuerdo a la función CC[1], a tiempos largos esta
	temperatura alcanza asintóticamente su nivel estacionario T2 y la temperatura en la barra tiende al estado
	estacionario, dado por

		f = 1 - x,	es decir,	T = T1 + (T2 - T1) * (1 - x)

	En referencia a la precisión del método numérico empleado, a medida que aumentamos el valor de N disminuye el valor de h (espaciado de la discretización espacial). A medida que h disminye debemos disminuir el valor del parámetro "pasork" (el espaciado de la discretización de la variable temportal) para que el sistema numérico esté bien definido y sea convergente.

	Aunque esto queda fuera del temario de FCI, y fuera de lo que se pide en la PEC, lo que sucede en este método numérico es que para que el sistema numérico se comporte bien a medida que h disminuye, en esta ecuación particular, debemos variar el parámetro "pasork" de forma proporcional a h^2.

	A medida que disminuimos los parámetros h y "pasork" (de manera coherente entre sí)  generamos una aproximación numérica cada vez más precisa. pero con un coste computacional (tiempo de cálculo) cada vez mayor. Esto es lo que sucede siempre que se hacen cálculos numéricos, y nos obliga a buscar constantemente un compromiso entre tiempo de cálculo y precisión, lo cual no siempre es fácil.

	Lo más recomendable para explorar el comportamiento del sistema estudiado en esta PEC es partir de un N no muy grande (p. ej. N = 10), y un paso temporal de (p. ej.) pasork = 0.01, e ir incrementando N poco a poco (multiplicando por 2 cada vez) y disminuyendo el paso poco a poco (dividiendo por 4 cada vez).
*/