-
Notifications
You must be signed in to change notification settings - Fork 4
/
Copy pathmimetadap.m
78 lines (67 loc) · 1.99 KB
/
mimetadap.m
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
function [t,u,H,ERROR,rechazo] = mimetadap(t0, tfin, x0, l0,hmin, TOL, metodo, orden, f, par)
% MÉTODO ADAPTATIVO - con un método monopaso y estimación del error en dos
% pasos sucesivos
%
% DATOS DE ENTRADA
% t0 tiempo inicial
% tfin tiempo final T
% x0 valor inicial, vector columna
% l paso inicial
% hmin paso minimo (por ejemplo hmin=10^{-5})
% TOL tolerancia para el error local (por ejemplo 10^{-3})
% metodo método a utilizar
% orden orden del método
% f(t,x) funcion de la E.D.O.
% par puede ser una variable de entrada con parametros,se incorpora al final de los
% datos de entrada, y si no los hay se define como una variable vacia par = []
%
% DATOS DE SALIDA
% t vector fila de tiempos t(n)
% u tabla de valores de x(t(n))
% H un vector fila de pasos
% ERROR un vector fila de errores
% rechazo el numero de veces que se rechaza el calculo
hmax = (tfin-t0)/10; % Por ejemplo
fac = 0.9; % Factor de seguridad
facmax = 5; % Factor de contracción/dilatación
rechazo =0;
l=l0;
ERROR=[];
H=[];
%Inicializamos t
t(1) = t0;
%inicializar la variable u como una tabla de ceros (la columna n+1 representa
%la aproximación de x(tn)), añadiendo el valor inicial
u(:,1) = x0;
%Iteramos
i=1;
while t(i)<tfin
%Calculamos el nodo tentativo
tten= t(i)+l;
if tten>tfin
tten=tfin;
l=tfin-t(i);
end
%if tten > tfin
%break;
%end
%Calculamos el valor tentativo
[~,x]=metodo(t(i), tten, 1, u(:,i), f, par);
%Calculamos el valor con dos pasos para estimar el error
[~,xx]=metodo(t(i), tten, 2, u(:,i), f, par);
%Estimamos el error (resta en valor absoluto / l)
ERR= norm(x(:,2)-xx(:,3))/l;
ERROR=[ERROR, ERR];
if ERR <= TOL
t(i+1)=tten;
H(i)=l;
u(:,i+1)=x(:,2);
i = i+1;
elseif l < hmin
error('Paso demasiado pequeño')
else
rechazo=rechazo+1;
end
l=min(hmax,l*min(facmax,fac*((TOL/ERR)^(1/orden))));
end
end