Metodos numericos

1.

(a) Cuando una partícula se mueve en una región donde hay un campo eléctrico y uno magnético, la fuerza total es la suma de la fuerza eléctrica y la magnética:

Que es la denominada fuerza de Lorentz. En el problema:
Entonces:

De la Segunda Ley de Newton, a:

O, equivalentemente:

(b) Considerando las condiciones iniciales:

Elegimos las variables de estado:
-Variables deestado:


-Ecuación de estado:
PVI1
(c) En laboratorio se implementó ya el método GBS, que aquí reproducimos:

PVI_GBS_2010_2.m

function [tv,yv]=PVI_GBS_2010_2(fun_f,dom,alpha,h,Nsig)

%input
% fun_f: función pendiente
% dom: vector dominio, dom=[t0,tf]
% alpha: condición inicial, y(t0)=alpha
% h: tamaño de paso, h=(t0-tf)/n, n es el número de particiones
%Nsig: Presición, número de cifras significativas de la aproximación
%
%output
% tv: vector de nodos, tv=[t0,t1,t2,...,tn]
% yv: vector de valores aproximados, yv=[y0,y1,y2,...,yn]

Tol = 0.5*10^(-Nsig); % La manera ya conocida de definir la Tolerancia
t0 = dom(1); tf = dom(2); % Establecemos los nodos extremos como los extremos
% del dominio
n = round((tf-t0)/h); % El númerode particiones
neq = length(alpha); % El número de ecuaciones del sistema
tv = zeros(n+1,1); % Inicializamos las variables que almacenarán los
yv = zeros(n+1,neq); % nodos y las aproximaciones
t = t0; tv(1,1) = t; % Iniciamos en el primer nodo
y = alpha; yv(1,1:neq) = y; % Usamos las condiciones iniciales establecidas

% Algoritmo BGS %

% Cálculo de y_(j+1)[h, Nsig]
for j =2:(n+1) % Se toma desde 2 pues MATLAB inicia en y(1) y no y(0)
k = 1; OK = 1; % Usamos un contador k y un criterio de parada OK
q1 = 2; q2 = 3; % Primeros términos de la sucesión de Bulirsch-Stoer
while OK
hb = h/q1; %tamaño de paso h barra
tb = t; % nodo sobre el que se trabaja
y0 = y; %y0 barra, valor inicial del intervalo a particionar
yb = y0+ hb*fun_f(tb,y0); % Método de Euler para y1 (ó y(2))
y1 = yb; % Almacenamos el valor obtenido
tb = tb + hb; % Avanzamos al siguiente “sub-nodo”
for i=3:(q1+1) % Aplicamos punto medio para yk,
yb = y0 + 2*hb*fun_f(tb,y1); % k=2,…q
y0 = y1; % Idem, comentarios anteriores en Euler
y1 = yb;
tb = tb+ hb;
end
yb = 1/2*(y0 + yb + hb*fun_f(tb,yb)); % Corrección de yb_q.
% o Fórmula de Gragg
if isequal(k,1)
F1(1,:) = yb; % Esto para la primera iteración
else
F2(1,:) = yb;
F2 = fun_Richardson(F1,F2,q1,q2); %devuelve la fila 2 comleta, el último término viene a ser en elemento en la diagonal
ifnorm(F2(k,:)-F1(k-1,:),inf) <= Tol*norm(F2(k,:),inf) % Es la norma infinita. También se puede escribir (k,:) por (k,1:neq)
y_G = F2(k,:); % Si se cumple la condición de presición,
% almacenamos el valor encontrado
F2 = []; % Reinicializamos F2
OK = 0; % Criterio de parada del bucle
else
aux = q1; % Si no usamos el siguienteqk de la sucesión
q1 = q2; % Esta es la fórmula de recurrencia de
q2 = 2*aux; % la sucesión de Bulirsch-Stoer
end
F1 = F2; % Actualizamos para generar una nueva fila
end
k = k+1; % Aumentamos en uno el contador
end %while
y = y_G; yv(j,1:neq) = y; % Actualizamos los valores antes de iniciar la
t= t + h; tv(j,1) = t; % siguiente iteración
end %for
%end_GBS

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function F2 = fun_Richardson(F1,F2,q1,q2)

k = size(F1,1);
p = 2; %p es el orden del método, y la fórmula de Gragg es de orden 2
for j = 2:(k+1) % pues ya se estableció F2(1,:) = yb
n = p*(j-1); % De la fórmula que se obtuvo en clase
F2(j,:) = ( q2^n*F2(j-1,:)...