## La sucesión de Fibonacci 

f[1]:1$ f[2]:1$
f[n]:=f[n-1]+f[n-2];

f[30];

lista1:makelist(f[n],n,1,20);

aa:(1+sqrt(5))/2;
bb:(1-sqrt(5))/2;
ff(n):=1/sqrt(5)*(aa^n-bb^n);

ff(11);

ratsimp(%);

f[11];

makelist(ratsimp(ff(k))-f[k],k,1,20);

fpprec:30;

N:20$
l:bfloat(ff(N+1)/ff(N));
l1:bfloat(aa);
l-l1;

coef:makelist(bfloat(ff(n+1)/ff(n)),n,1,N)$

wxplot2d([[discrete, coef], l1],[x,1,N], [y,1,2],[style, points, lines],
          [point_type,diamond], [color,blue,red],[legend,false])$

N:20$
nn:makelist(n,n,1,N)$
coef:makelist(bfloat(ff(n+1)/ff(n)),n,1,N)$
er:makelist(bfloat(expand(ff(n+1))/expand(ff(n))-l1),n,1,N)$
transpose(matrix(nn,coef,er));

/* Solucion analitica: Paquete solve_rec */
/* Si quieres sabr que hace ejecuta la orden:  demo(solve_rec) */
load("solve_rec")$

ec:fib1[n]=fib1[n-1]+fib1[n-2];

solve_rec(ec,fib1[n]);

define(fiba[n],second(solve_rec(ec,fib1[n],fib1[1]=1,fib1[2]=1)));

fiba[n];

fiba[3];radcan(%);

makelist(radcan(fiba[n]),n,1,20)-lista1;

fiba[n+1]/fiba[n]$
radcan(%);

/* Este resultado es sospechoso: no te fies nunca de un ordenador */
limit(%,n,inf);

define(cociente[n],radcan(fiba[n+1]/fiba[n]));

factor(%);

float(cociente[100]);

float((1+sqrt(5))/2);

makelist([n,float(cociente[n])],n,1,15);

wxplot2d([discrete,%], [style, points], [point_type, bullet]);

define(sus(n),cociente[n]);

float(cociente[100.0]);

limit(sus(n),n,+inf);

/* Intentando aclarar que está ocurriendo */ 

assume(a>b);

define(eje(n),a^n+b^n);

ratsimp(eje(n+1)/eje(n));

limit(%,n,inf);

define(eje(n),(-a)^n+(-b)^n);

ratsimp(eje(n+1)/eje(n));

limit(%,n,inf);

kill(all);

## El método de Newton para calcular raíces 1669 publicado en 1711

### Ejemplo 1: Ceros de un polinomio de grado 3

kill(p,dp,a)$
p(x):=x*(x-1)*(x-2);

wxplot2d(p(x),[x,-1/2,5/2]);

define(dp(x),diff(p(x),x));

a[1]:4$ 
a[n]:=a[n-1]-p(a[n-1])/dp(a[n-1]);

a[2];

N:10$
lista:makelist(float(a[n]),n,1,N);

wxplot2d([[discrete, lista], 2],[x,1,N], [y,1.5,4],[style, points, lines],
          [point_type,diamond], [color,blue,red],[legend,false])$

wxplot2d([[discrete, lista], 2],[x,1,N], [y,1.5,4],[style, points, lines])$

nn:makelist(n,n,1,N)$
aa:makelist(float(a[n]),n,1,N)$
er:makelist(float(2-a[n]),n,1,N)$
transpose(matrix(nn,aa,er));

### Ejemplo 2. El sen(x)

kill(p,dp,a)$
p(x):=sin(x);
define(dp(x),diff(p(x),x));
a[1]:2$ 
a[n]:=a[n-1]-p(a[n-1])/dp(a[n-1]);

N:10$
lista:makelist(float(a[n]),n,1,N);

wxplot2d([[discrete, lista], π],[x,1,N], [y,2,4],[style, points, lines],
          [point_type,diamond], [color,blue,red],[legend,false])$

nn:makelist(n,n,1,N)$
aa:makelist(float(a[n]),n,1,N)$
er:makelist(float(π-a[n]),n,1,N)$
transpose(matrix(nn,aa,er));

### Ejemplo 3: Ceros de un polinomio de grado 5

kill(p,dp,a)$
p(x):=x^5+3*x^2-1;

wxplot2d(p(x),[x,-2,2],[y,-10,10]);

define(dp(x),diff(p(x),x));
a[1]:1.0$ 
a[n]:=a[n-1]-p(a[n-1])/dp(a[n-1]);

a[30];

N:10$
lista:makelist(float(a[n]),n,1,N);

wxplot2d([[discrete, lista]],[x,1,N], [style, points],
          [point_type,diamond], [color,blue,red],[legend,false])$

nn:makelist(n,n,2,N)$
aa:makelist(float(a[n]),n,2,N)$
er:makelist(float(a [n]-a[n-1]),n,2,N)$
transpose(matrix(nn,aa,er));

kill(a)$
a[1]:-0.2$ /* probar a[1]:0 y 0.1 */
a[n]:=a[n-1]-p(a[n-1])/dp(a[n-1]);
N:10$
lista:makelist(float(a[n]),n,1,N);
nn:makelist(n,n,2,N)$
aa:makelist(float(a[n]),n,2,N)$
er:makelist(float(a [n]-a[n-1]),n,2,N)$
transpose(matrix(nn,aa,er));

kill(a)$
a[1]:-2.0$ /* probar a[1]:-1.0 */
a[n]:=a[n-1]-p(a[n-1])/dp(a[n-1]);
N:10$
lista:makelist(float(a[n]),n,1,N);
nn:makelist(n,n,2,N)$
aa:makelist(float(a[n]),n,2,N)$
er:makelist(float(a [n]-a[n-1]),n,2,N)$
transpose(matrix(nn,aa,er));

## Modelos discretos de poblaciones

/* Modelos discretos: ecuaciones en diferencias */

kill(all)$
load(draw);

### El modelo de Malthus discreto 

/* Considerar r>0 y r<0 */
kill(p)$
p[1]:20; r:.1;
p[n]:=(r+1)*p[n-1];

p[5];
pobm1:makelist(p[k],k,1,100)$

wxdraw2d(point_type = 5,point_size = 0.5, color="#0000ff", points(pobm1))$

### El modelo de Verhulst  p[n]:=r*p[n-1]/(K+p[n-1])

kill(all)$
solve(p='r*p/('K+p),p);    /*  hacer para p[1]<K-r y p[1]>K-r */

/* Considerar r>K y r<K */
kill(all)$ /* r>K */
load ("descriptive")$
p[1]:8; r:20.0; K:10;  /*  hacer para p[1]<K-r */
p[n]:=r*p[n-1]/(K+p[n-1]);

p[5];
pobm2:makelist(p[k],k,1,10);
makelist(p[k+1]/p[k],k,1,10);
wxplot2d([p,r*p/(K+p),r-K],[p,0,30]);      

mx:smax(pobm2)$ mi:smin(pobm2)$
wxdraw2d(point_type = 7,point_size = 1, color="#0000ff", points(pobm2),yrange=[mi-1,mx+1],explicit(r-K,x,0,10))$

kill(all)$ /* r>K */          
load ("descriptive")$
p[1]:30; r:20.0; K:10;  /*  hacer para p[1]>K-r */
p[n]:=r*p[n-1]/(K+p[n-1]);

p[5];
pobm2:makelist(p[k],k,1,10);
makelist(p[k+1]/p[k],k,1,10);
wxplot2d([p,r*p/(K+p),r-K],[p,0,30]);

mx:smax(pobm2)$ mi:smin(pobm2)$
wxdraw2d(point_type = 7,point_size = 1, color="#0000ff", points(pobm2),yrange=[mi-1,mx+1],explicit(r-K,x,0,10))$

kill(p)$  /*  r<K  */
p[1]:1; r:2.0; K:10;
p[n]:=r*p[n-1]/(K+p[n-1]);
p[5];
pobm2:makelist(p[k],k,1,10);

makelist(p[k+1]/p[k],k,1,10);
wxplot2d([p,r*p/(K+p),0,r-K],[p,0,30]);

mx:smax(pobm2)$ mi:smin(pobm2)$
wxdraw2d(point_type = 7,point_size = 1, color="#0000ff", points(pobm2),yrange=[mi-1,mx+1],explicit(0,x,0,10))$

limit('r*p/('K+p),p,inf);

### Sucesiones por recurrencia: método general

Estudiando la relación de recurrencia x(n+1)=x(n)(4/3-x(n))

kill(all);

wxplot2d([x*(4/3-x),x],[x,-1,2]);

wxplot2d([x*(4/3-x),x],[x,0,0.4],[legend,false]);

der:diff(x*(4/3-x),x);wxplot2d(%,[x,-1,2]);

solve(der,x);ev(x*(4/3-x),%);

float(4/9);

float(4/3);

kill(x)$
x[1]:0.1;
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,4);
lis:makelist(x[k],k,1,30)$
wxplot2d([discrete,lis]);

kill(x)$
x[1]:2/3+.1; /* probar 2/3-0.1 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,4);
lis:makelist(x[k],k,1,30)$
wxplot2d([discrete,lis]);

kill(x)$
x[1]:-1.1; /* probar -1.0 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,10);
wxplot2d([discrete,%],[y,-10000,1]);

kill(x)$
x[1]:1; /* probar 2 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,10);
wxplot2d([discrete,%]);

kill(x)$
x[1]:1.2; /* probar 2 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,20);
wxplot2d([discrete,%]);

wxplot2d([x,x*(4/3-x)],[x,-1/2,5/3]);

wxplot2d([x,x*(4/3-x)],[x,-.02,1/2.9]);

solve(x-x*(4/3-x)=0,x);

Estudiando el modelo de Verhulst

kill(all)$
r*p[n]/(K+p[n])-r*p[n-1]/(K+p[n-1]);

ratsimp(%);

### El modelo logístico discreto 

load(draw);

kill(r,ld)$
r:2;
ld[1]:0.4;
ld[n]:=r*ld[n-1]*(1 - ld[n - 1]);

pl:makelist(ld[n],n,1,300)$

wxplot2d([discrete,pl] ,[style, points],[y,0.35,0.6],
          [point_type,diamond], [color,blue,red],[legend,false]);

wxdraw2d(yrange=[0.35,0.6],point_type = 7,point_size = 1, color="#0000ff", points(pl));

kill(r,ld,pl)$
r:3.01;
ld[1]:0.4;
ld[n]:=r*ld[n-1]*(1 - ld[n - 1]);
ld[10];
pl:makelist(ld[n],n,1,300)$

wxplot2d([discrete,pl] ,[style, points],
          [point_type,diamond], [color,blue,red],[legend,false]);

wxdraw2d(yrange=[0.35,0.8],point_type = 7,point_size = 1, color="#0000ff",points(pl))$

3.46;

wxdraw2d(point_type = 5,point_size = 0.5, color="#0000ff", points(pl))$

wxdraw2d(yrange=[0.35,0.9],point_type = 7,point_size = 1, color="#0000ff",points(pl))$

kill(r,ld,pob1)$
r:3.56;
ld[1]:0.4;
ld[n]:=r*ld[n-1]*(1 - ld[n - 1]);
ld[10];
pob1:makelist(ld[n],n,1,1000)$
wxdraw2d(point_type = 7,point_size = 0.5, color="#0000ff", points(pob1),yrange=[0.25,.95])$

wxdraw2d(yrange=[0.3,1],point_type = 8,point_size =0.4, color="#0000ff",points(pob1))$

kill(r,ld,pob1)$
r:3.566;
ld[1]:0.4;
ld[n]:=r*ld[n-1]*(1 - ld[n - 1]);
ld[10];
pob1:makelist(ld[n],n,1,2000)$
wxdraw2d(point_type = 7,point_size = 0.5, color="#0000ff", points(pob1),yrange=[0.3,.95])$

wxdraw2d(yrange=[0.3,1.],point_type = 7,point_size =0.5, color="#0000ff",points(pob1))$

/* OPCIONAL: Otra forma de hacerlo */
load("dynamics")$

/*  r = 2 */
evolution(2*y*(1-y), 0.4, 100);

/*  r = 3.46 */
evolution(3.46*y*(1-y), 0.4, 300)$

/*  r = 3.46 */
evolution(3.566*y*(1-y), .4, 1000)$;

/*  r = 3.6 */
evolution(3.6*y*(1-y), 0.4, 1000)$

kill(x,y);
orbits(x*y-x*y^2, 0.4, 50, 100, [x, 2.6, 3.9, 0.01], [style, dots]);

## Sucesiones por recurrencia: método general

Estudiando la relación de recurrencia x(n+1)=x(n)(4/3-x(n))

kill(all);

wxplot2d([x*(4/3-x),x],[x,-1,2]);

wxplot2d([x*(4/3-x),x],[x,0,0.4],[legend,false]);

der:diff(x*(4/3-x),x);wxplot2d(%,[x,-1,2]);

solve(der,x);ev(x*(4/3-x),%);

float(4/9);

float(4/3);

kill(x)$
x[1]:0.1;
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,4);
lis:makelist(x[k],k,1,30)$
wxplot2d([discrete,lis]);

kill(x)$
x[1]:2/3+.1; /* probar 2/3-0.1 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,4);
lis:makelist(x[k],k,1,30)$
wxplot2d([discrete,lis]);

kill(x)$
x[1]:-1.1; /* probar -1.0 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,10);
wxplot2d([discrete,%],[y,-10000,1]);

kill(x)$
x[1]:1; /* probar 2 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,10);
wxplot2d([discrete,%]);

kill(x)$
x[1]:1.2; /* probar 2 */
x[n]:=x[n-1]*(4/3-x[n-1]);
makelist(x[k],k,1,20);
wxplot2d([discrete,%]);

wxplot2d([x,x*(4/3-x)],[x,-1/2,5/3]);

wxplot2d([x,x*(4/3-x)],[x,-.02,1/2.9]);

solve(x-x*(4/3-x)=0,x);