# Tema 12
12-4  Introducción a los métodos numéricos
          de resolución de ecuaciones diferenciales 
          ordinarias y sistemas

Ejemplo 12.4.1

kill(all)$

EDO: t*(t-1)*diff(y(t), t) - (1-2*t)*y(t) + y(t)^2 = 2*t;

ode2(EDO, y(t), t);

EDO_1b : (t-1)*t*Dy(t)+y(t)^2-(1-2*t)*y(t)=2*t;

solve( [EDO_1b]  ,  [Dy(t)] );

fpprintprec : 5$

soluc1: rk( [-(y^2+(2*t-1)*y-2*t)/(t^2-t) ], [y], [2] , [t, 2, 5, 0.1] );

wxplot2d( [ [discrete, soluc1] ] , [style, [points, 1, 1, 2] ]  );

wxplot2d( [ [discrete, soluc1] ] , [style, [lines, 2, 1 ] ]  );

Ejemplo 12.4.2

kill(all)$

EDO: (100 - t^6)*diff(y(t), t) - cos(t) = 0;

ode2(EDO, y(t), t);

fpprintprec : 5 $

soluc1: rk( [ cos(t)/(100 - t^6) ], [y], [0.01] , [t, 0, 1,  0.05] );

wxplot2d( [ [discrete, soluc1] ] , [style, [points, 1, 1, 2] ]  );

Ejemplo 12.4.3

kill(all)$
load("dynamics")$   numer:true$     grind:true$

edo1:  'diff(x1, t)=   -0.2*x1 +  0.3*x2;
edo2:  'diff(x2, t)=   0.1*x1 - 0.4*x2;

x1_0:1;          x2_0:3;
t_ini=t0:0.0;     t_fin=t1:10.0;        Pas=p:0.5;     n_int=n_int:(t1 - t0)/0.5;

fpprintprec : 5 $

solucio_LV: rk([ - 0.2*x1 + 0.3*x2 , 0.1*x1 - 0.4*x2], [x1, x2], [x1_0, x2_0], [t, t0, t1, p]);

M1:matrix(
[0.0,1,3],[0.5,1.2961,2.5088],[1.0,1.5014,2.1179],[1.5,1.6374,1.8055],[2.0,1.7204,1.5545],
[2.5,1.7634,1.3518],[3.0,1.7762,1.1871],[3.5,1.7665,1.0523],[4.0,1.7403,0.941],
[4.5,1.7021,0.848],[5.0,1.6554,0.771],[5.5,1.603,0.705],[6.0,1.5468,0.648],[6.5,1.4886,0.6],
[7.0,1.4294,0.557],[7.5,1.3701,0.519],[8.0,1.3113,0.486],[8.5,1.2537,0.456],
[9.0,1.1975,0.429],[9.5,1.1429,0.404],[10.0,1.0902,0.381]
)$

M1T:transpose(M1)$
T:first(M1T)$       Y1:second(M1T)$      Y2:third(M1T)$

List1:makelist([T[n],Y1[n]],  n, 1, 1 + n_int)$
List2:makelist([T[n],Y2[n]],  n, 1, 1 + n_int)$

wxplot2d( 
[ [discrete, List1] , [discrete, List2] ] , 
[style, [lines, 2, 1] , [lines, 2, 2] , [lines, 2, 3] ] 
);

Ejemplo 12.4.4

kill(all)$
load("dynamics")$   numer:true$     grind:true$

edo1:'diff(N1, t)=   a1*N1 -  a2*N1*N2;
edo2:'diff(N2, t)= - b1*N2 + b2*N1*N2;

N1_0:30;          N2_0:6;
a1:0.5654;        a2:0.0245;        b1:0.8267;        b2:0.0219;

t_ini=t0:0.0;     t_fin=t1:20.0;        Pas=p:0.5;     n_int=n_int:(t1 - t0)/0.5;

fpprintprec : 5 $

solucio_LV: rk([ a1*N1 -  a2*N1*N2, - b1*N2 + b2*N1*N2], [N1, N2], [N1_0, N2_0], [t, t0, t1, p]);

M1:matrix(
[0.0,30,6],[0.5,37.062,5.721],[1.0,45.805,5.9473],[1.5,56.232,6.8679],[2.0,67.82,8.9548],
[2.5,78.782,13.243],[3.0,84.883,21.613],[3.5,79.734,35.646],[4.0,61.802,51.611],[4.5,40.965,59.755],
[5.0,26.386,56.753],[5.5,18.405,47.734],[6.0,14.477,37.697],[6.5,12.794,28.905],[7.0,12.455,21.932],
[7.5,13.065,16.669],[8.0,14.488,12.811],[8.5,16.727,10.046],[9.0,19.871,8.1114],[9.5,24.073,6.8174],
[10.0,29.533,6.0399],[10.5,36.474,5.7259],[11.0,45.087,5.9101],[11.5,55.396,6.7669],[12.0,66.937,8.7396],
[12.5,78.055,12.808],[13.0,84.734,20.793],[13.5,80.615,34.424],[14.0,63.416,50.564],[14.5,42.366,59.56],
[15.0,27.223,57.265],[15.5,18.835,48.491],[16.0,14.678,38.423],[16.5,12.867,29.504],[17.0,12.445,22.393],
[17.5,12.991,17.011],[18.0,14.355,13.06],[18.5,16.529,10.222],[19.0,19.601,8.2328],[19.5,23.717,6.8958],
[20.0,29.075,6.0823]
)$

M1T:transpose(M1)$
T:first(M1T)$       Y1:second(M1T)$      Y2:third(M1T)$     M:Y1+Y2$

List1:makelist([T[n],Y1[n]],  n, 1, 1 + n_int)$
List2:makelist([T[n],Y2[n]],  n, 1, 1 + n_int)$
List3:makelist([T[n],M[n]],  n, 1, 1 + n_int)$

wxplot2d( 
[ [discrete, List1] , [discrete, List2] , [discrete, List3] ] , 
[style, [lines, 2, 1] , [lines, 2, 2] , [lines, 2, 3] ] 
);