forked from mx-psi/mn2
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathMétodosT3.mac
151 lines (131 loc) · 5.1 KB
/
MétodosT3.mac
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
/****************************************************************\
** Métodos requeridos en los ejercicios del Tema 3 **
** Autores: Pablo Baeyens Fernández y José Manuel Muñoz Fuentes **
\****************************************************************/
/* Métodos de Taylor de orden r para resolver el PVI */
/* Incluye el método de Euler, que es el de Taylor de orden 1 */
pvi_taylor(a, b, n, y0, f, r) := block([t,h,F,T,u],
u[0]: y0,
t: a,
h: (b-a)/n,
F[0](x,y) := f(x,y),
F[n](x,y) := diff(F[n-1](x,y),x)+diff(F[n-1](x,y),y)*f(x,y),
define(T(x,y), f(x,y) + sum(h^k*F[k](x,y)/(k+1)!, k, 1, r-1)),
for i:1 thru n do (
u[i]: float(u[i-1] + h*T(t, u[i-1])),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$
pvi_euler(a, b, n, y0, f) := pvi_taylor(a, b, n, y0, f, 1)$
/* Método de Euler mejorado para resolver el PVI */
/* Es un método de Runge-Kutta de orden 2 */
pvi_euler_mejorado(a, b, n, y0, f) := block([u,h,t],
u[0]: y0,
t: a,
h: (b-a)/n,
for i:1 thru n do (
u[i]: float(u[i-1] + h*f(t+h/2,u[i-1]+ h/2*f(t,u[i-1]))),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$
/* Método del punto medio para resolver el PVI. Nótese que se trata de un método de dos pasos: */
/* u1 es una aproximación del valor en el segundo nodo obtenida con otro método para resolver el PVI */
pvi_punto_medio(a, b, n, u0, u1, f) := block([u,t,h],
u[0]: u0,
u[1]: u1,
h: (b-a)/n,
t: a+h,
for i:2 thru n do(
u[i]: float(u[i-2] + 2*h*f(t, u[i-1])),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$
/* Método del trapecio o de Crank-Nicholson para resolver el PVI */
/* Es implícito: es necesario encontrar el cero de una función */
/* Esto lo hacemos con la función newton, que aplica el método de Newton-Raphson */
/* Tomamos como aproximación inicial el valor en el nodo anterior */
load(newton1)$
pvi_crank_nicholson(a, b, n, u0, f) := block([u,t,h],
u[0]: u0,
h: (b-a)/n,
t: a+h,
for i:1 thru n do(
u[i]: float(newton(x-(u[i-1] + h/2*(f(t-h, u[i-1] + f(t, x)))), x, u[i-1], 10^(-4))),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$
/* Método de Heun para resolver el PVI */
/* Es un método de Runge-Kutta de orden 2 */
pvi_heun(a, b, n, y0, f) := block([u,h,t,K],
u[0]: y0,
t: a,
h: (b-a)/n,
for i:1 thru n do (
K: f(t, u[i-1]),
u[i]: float(u[i-1] + h/2*(K + f(t+h, u[i-1] + h*K))),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$
/* El método de Runge-Kutta de orden 4 */
/* Es un método de Runge-Kutta de orden 4 aparentemente muy usado */
pvi_runge_kutta_4(a, b, n, y0, f) := block([u,h,t,K1,K2,K3],
u[0]: y0,
t: a,
h: (b-a)/n,
for i:1 thru n do (
K1: f(t, u[i-1]),
K2: f(t+h/2, u[i-1] + h/2*K1),
K3: f(t+h/2, u[i-1] + h/2*K2),
K4: f(t+h, u[i-1] + h*K3),
u[i]: float(u[i-1] + h/6*(K1 + 2*K2 + 2*K3 + K4)),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$
/* Método de Adams-Bashforth. Se trata de un método multipaso que estima el valor en el */
/* siguiente nodo usando la integral del polinomio que interpola los k valores anteriores */
/* El método de k pasos se calcula una sola vez para cada k y cada ejecución del programa */
/* Requiere pasar como parámetro un método para obtener los valores en los primeros nodos */
/* Usando formula_adams_bashforth[k] se puede visualizar la fórmula particular para k pasos */
formula_adams_bashforth[k](h, ftu, t) := factor(integrate(pn(ftu, makelist(t-(k-i)*h, i, 0, k-1)), x, t-h, t))$
/* Para obtener el polinomio de interpolación utilizamos los métodos del tema anterior */
load(sconcat(pathname_directory(load_pathname), "MétodosT2.mac"))$
pvi_adams_bashforth(a, b, n, y0, metodo_inicial, f, k) := block([u,h,t,m],
h: (b-a)/n,
t: a+h*k,
m: metodo_inicial(a, t, k, y0, f),
f_(t) := f(t, u[(t-a)/h]),
for i:0 thru k-1 do u[i]: m[i+1][2],
for i:k thru n do (
u[i]: float(u[i-1] + formula_adams_bashforth[k](h, f_, t)),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$
/* Método de Adams-Moulton. Se trata de un método multipaso implícito que estima el valor */
/* en el siguiente nodo con la integral del polinomio que interpola los k valores previos */
/* El método de k pasos se calcula una sola vez para cada k y cada ejecución del programa */
/* Requiere pasar como parámetro un método para obtener los valores en los primeros nodos */
/* Es implícito: requiere encontrar el cero de una función */
/* Esto lo hacemos con la función newton, que aplica el método de Newton-Raphson */
/* Tomamos como aproximación inicial el valor en el nodo anterior */
/* Usando formula_adams_moulton[k] se puede visualizar la fórmula particular para k pasos */
formula_adams_moulton[k](h, ftu, t) := factor(integrate(pn(ftu, makelist(t-(k-i)*h, i, 0, k)), x, t-h, t))$
pvi_adams_moulton(a, b, n, y0, metodo_inicial, f, k) := block([u,h,t,m],
h: (b-a)/n,
t: a+h*k,
m: metodo_inicial(a, t, k, y0, f),
f_(t) := f(t, u[(t-a)/h]),
for i:0 thru k-1 do u[i]: m[i+1][2],
for i:k thru n do (
u[i]: x,
u[i]: float(newton(x-(u[i-1] + formula_adams_moulton[k](h, f_, t)), x, u[i-1], 10^(-4))),
t: t+h
),
makelist([a+h*i, u[i]], i, 0, n)
)$