Zeros de Funções
Método da Bissecção:
Os arcos foram adicionados para destacar o intervalo de cada iteração.
Método da Posição Falsa:
Os arcos foram adicionados para destacar o intervalo de cada iteração.
Método do Ponto Fixo:
Método de Newton:
Método da Secante:
Sistemas Lineares - Métodos Diretos
Sistema Triangular Inferior (substituição direta/progressiva)
function X=subdireta(A,b) n=length(A); x(1)=b(1)/A(1,1) for i=2:n s=0; for j=1:i-1 s=s+A(i,j)*x(j); end x(i)=(b(i)-s)/A(i,i); end X=x; end
Sistema Triangular Superior (substituição regressiva)
function X=subreg(A,b) x(n)=b(n)/A(n,n) for i=n-1:-1:1 s=0; for j=i+1:n s=s+A(i,j)*x(j); end x(i)=(b(i)-s)/A(i,i); end X=x; end
Eliminação de Gauss
function [U c]=Gauss(A,b) n=length(A); for j=1:n-1 for i=j+1:n m(i,j)=A(i,j)/A(j,j); b(i)=b(i)-m(i,j)*b(j); for k=1:n A(i,k)=A(i,k)-m(i,j)*A(j,k); end end end U=A; c=b; end
Eliminação de Gauss com pivoteamento Parcial
function [U c]=GaussP(A,b) n=length(A); for j=1:n-1 w=abs(A(j,j)) for i=j:n if abs(A(i,j))>w w=abs(A(i,j)); r=i; end end v=A(j,:); A(j,:)=A(r,:); A(r,:)=v; for i=j+1:n m(i,j)=A(i,j)/A(j,j); b(i)=b(i)-m(i,j)*b(j); for k=1:n A(i,k)=A(i,k)-m(i,j)*A(j,k); end end end U=A; c=b; end
Sistemas Lineares - Métodos Iterativos
Métodos de Jacobi (preto) e Gauss-Seidel (vermelho):Método de Jacobi
function X=jacobi(A,b,x,kmax,tol) k=0; dr=tol+1; while k<=kmax && dr>tol k=k+1; for i=1:n s=0; for j=1:n s=s+A(i,j)*x(j); end s=s-A(i,i)*x(i); xx(i)=(b(i)-s)/A(i,i); end dr=norm(xx-x,inf); % maximo dos módulos das diferenças das componentes de xx e x x=xx; end X=x; end
Método de Jacobi (versão matricial)
function X=Jacobi(A,b,x,kmax,tol) k=0; dr=tol+1; D=diag(A); % matriz diagonal dos elementos da diagonal principal de A Di=inv(D); % inversa de D Ad=A-D; while k<=kmax && dr>tol k=k+1; xx=Di*(b-Ad*x); dr=norm(xx-x,inf); x=xx; end X=x; end
Método de Gauss-Seidel
function X=gseidel(A,b,x,kmax,tol) k=0; dr=tol+1; while k<=kmax && dr>tol k=k+1; s=0; for j=2:n s=s+A(1,j)*x(j); end xx(1)=(b(1)-s)/A(1,1); for i=2:n s=0; for j=1:i-1 s=s+A(i,j)*xx(j); end if i<=n-1 for j=i+1:n s=s+A(i,j)*x(j); end end xx(i)=(b(i)-s)/A(i,i); end dr=norm(xx-x,inf); % maximo dos módulos das diferenças das componentes de xx e x x=xx; end X=x; end
Método de Gauss-Seidel (versão matricial)
function X=GSeidel(A,b,x,kmax,tol) k=0; dr=tol+1; L=tril(A); % matriz triangular inferior com L(i,j)=A(i,j) se i>=j e L(i,j)=0, caso contrário U=A-L; I=eye(size(A)); % matriz identidade % Cálculo da Inversa de L y=subdireta(L,I(:,1)); Li=y; for k=2:n y=subdireta(L,I(:,k)); Li=[Li y]; end % while k<=kmax && dr>tol k=k+1; xx=Li*(b-U*x); dr=norm(xx-x,inf); x=xx; end X=x; end
Critério das Linhas
function Alpha=criteriolinhas(A) for i=1:n s=0 for j=1:n s=s+abs(A(i,j)); end alpha(i)=(s-abs(A(i,i))/abs(A(i,i)); end Alpha=max(alpha); end
Critério de Sassenfeld
function Beta=sassenfeld(A) s=0; for j=2:n s=s+abs(A(1,j)); end beta(1)=s/abs(A(1,1)); for i=2:n s=0 for j=1:i-1 s=s+beta(j)*abs(A(i,j)); end if i<=n-1 for j=i+1:n s=s+abs(A(i,j)); end end beta(i)=s/abs(A(i,i)); end Beta=max(beta); end
Sistemas não-Lineares
Métodos de Newton (preto) e Newton Modificado (verde)
Método de Newton
function X=Newton(F,J,x,tol1,tol2,kmax) k=0; Dr=tol1+1; Fx=tol2+1; while k<=kmax && Dr>tol1 && Fx>tol2 k=k+1; s=GaussP(J(x),-F(x)); x=x+s; Dr=norm(s,inf); Fx=norm(F(x),inf); end X=x; end
Método de Newton Modificado
function X=Newtonmod(F,J,x,tol1,tol2,kmax) k=0; Dr=tol1+1; Fx=tol2+1; [L U P]=lu(J(x)); while k<=kmax && Dr>tol1 && Fx>tol2 k=k+1; y=subdireta(L,-P*F(x)); s=subreg(U,y); x=x+s; Dr=norm(s,inf); Fx=norm(F(x),inf); end X=x; end
Métodos de Broyden
function X=Broyden(F,B,x,tol1,tol2,kmax) k=0; Dr=tol1+1; Fx=F(x); while k<=kmax && Dr>tol1 && norm(Fx,inf)>tol2 k=k+1; s=GaussP(B,-Fx); x=x+s; Fx=F(x); B=B+Fx*s'/(s'*s); Dr=norm(s,inf); end X=x; end
Interpolação Polinomial
Forma de Newton: diferenças divididas e intervalos encaixantes
function fx=newtoninterpolation(x,y,t) n=length(x); %Calculo da matriz de diferencas divididas for i=1:n C(i,1)=y(i); end for j=2:n for i=1:n+1-j C(i,j)=(C(i+1,j-1)-C(i,j-1))/(x(i+j-1)-x(i)); end end %Calculo do polinomio em t na forma de parenteses encaixados d=C(1,1:n);%=C(1,:) fx=0; for k=n:-1:1 fx=d(k)+(t-x(k))*fx; end end
Exemplo de interpolação de uma função a partir de valores tabelados
x=0:0.2:2; y=cos(pi*x); tt=-1:0.01:2; for i=1:length(tt) p(i)=newtoninterpolation(x,y,tt(i)); end plot(tt,p,'r-','LineWidth',2) hold on z=cos(pi.*tt); plot(tt,z,'k--','LineWidth',2) plot(x,y,'bo','LineWidth',2) hold off grid on legend('interpolação','função','nós') erro=norm(p-z,inf) %calculado sobre os pontos tt(i)
Fenômeno de Runge e Nós de Chebyshev
Representação do fenômeno de Runge para pontos igualmente espaçados (laranja) e interpolação nos nós de Chebyshev para a função \[f(x)=\dfrac{1}{1+25x^2}\] Em verde, polinômio interpolador com nós: \[x_k=\dfrac{x_0+x_n}{2}-\dfrac{x_n-x_0}{2}\cos\left(\dfrac{k\pi}{n}\right)\] Em azul, polinômio interpolador com nós: \[x_k=\dfrac{x_0+x_n}{2}-\dfrac{x_n-x_0}{2}\cos\left(\dfrac{(2k+1)\pi}{2(n+1)}\right)\]
Interpolação por partes
Polinômio interpolador por partes quadrático (verde), polinômio interpolador de grau 6 (laranja) e spline cúbica (azul)
Polinômio interpolador por partes linear e polinômio interpolador de grau n (laranja)
Spline Cúbica x interpolação polinomial x linear por partes
Método dos Quadrados Mínimos
Código Matlab:
x=[-1;-0.75;-0.6;-0.5;-0.3;0;0.2;0.4;0.5;0.7;1]; y=[2.05; 1.15; 0.45; 0.4; 0.5; 0; 0.2; 0.6; 0.51;1.2;2.05]; g1=@(x) x.^2; g2=@(x) x; g3=@(x) 1; G=[g1(x) g2(x) ones(length(x))]; A=G'*G; b=G'*y; alpha=A\b; phi=@(x,c) c(1)*g1(x)+c(2)+c(3)*g3(x); J=@(c) sum((phi(x(:),c)-y(:)).^2); Jmin=J(alpha) %gráficos c=A\b; t=-1.3:0.01:1.3; plot(t,phi(t,alpha),'b-','LineWidth',3) hold on plot(x,y,'ro','LineWidth',5) grid on
Integração Numérica:
Quadratura Gaussiana
a=0;
b=1;
n=3;
f=@(x) exp(x);
[x, w] = Gauss(n);
x=0.5*(b+a+x*(b-a));
I=f(x)'*w;
I=0.5*(b-a)*I; %resposta final
function [x, w] = Gauss(n)
% Generates the abscissa and weights for a Gauss-Legendre quadrature.
% Reference: Numerical Recipes in Fortran 77, Cornell press.
x = zeros(n,1); % Preallocations.
w = x;
m = (n+1)/2;
for ii=1:m
z = cos(pi*(ii-.25)/(n+.5)); % Initial estimate.
z1 = z+1;
while abs(z-z1)>eps
p1 = 1;
p2 = 0;
for jj = 1:n
p3 = p2;
p2 = p1;
p1 = ((2*jj-1)*z*p2-(jj-1)*p3)/jj; % The Legendre polynomial.
end
pp = n*(z*p1-p2)/(z^2-1); % The L.P. derivative.
z1 = z;
z = z1-p1/pp;
end
x(ii) = -z; % Build up the abscissas.
x(n+1-ii) = z;
w(ii) = 2/((1-z^2)*(pp^2)); % Build up the weights.
w(n+1-ii) = w(ii);
end
end
Polinômios de Legendre (Octave)
a=0; b=1; n=3; pkg load miscellaneous c=legendrepoly(n); x=roots(c); x=0.5*(b+a+x*(b-a)); % calcular w
Problemas de Valor Inicial
\[y'(t)=f(t,y(t))\;,\quad y(t_0)=y_0\]
f=@(t,y) 3*y*(1-y)-y.^2/(1+y.^2);
t0=0; tf=10;
y(1)=0.1;
h=0.25;
n=(tf-t0)/h;
t=t0:h:tf;
for i=1:n
y(i+1)=y(i)+h*f(t(i),y(i));
end
plot(t,y,'-*','LineWidth',2)
legend(['h=',num2str(h)])
ylim([0 1.2])
grid on
title('Euler')
Método de Série de Taylor de ordem 2
f=@(t,y) 3*y*(1-y)-y.^2/(1+y.^2);
%derivadas
ft=@(t,y) 0;
fy=@(t,y) 3-6*y-2*y./(1+y.^2).^2;
t0=0; tf=10;
y(1)=0.1;
h=0.25;
n=(tf-t0)/h;
t=t0:h:tf;
for i=1:n
y(i+1)=y(i)+h*f(t(i),y(i))+((h^2)/2)*(ft(t(i),y(i))+fy(t(i),y(i))*f(t(i),y(i)));
end
plot(t,y,'-*','LineWidth',2)
legend(['h=',num2str(h)])
ylim([0 1.2])
grid on
title('Taylor ordem 2')
Método de Heun (Runge-Kutta de 2ª ordem)
f=@(t,y) 3*y*(1-y)-y.^2/(1+y.^2);
%derivadas
t0=0; tf=10;
y(1)=0.1;
h=0.25;
n=(tf-t0)/h;
t=t0:h:tf;
for i=1:n
K1=f(t(i),y(i));
K2=f(t(i)+h,y(i)+h*K1);
y(i+1)=y(i)+h*(K1+K2)/2;
end
plot(t,y,'-*','LineWidth',2)
grid on
title('Método de Heun (Runge-Kutta de ordem 2)')
Método de Runge-Kutta de 4ª ordem
f=@(t,y) 3*y*(1-y)-y.^2/(1+y.^2);
%derivadas
t0=0; tf=10;
y(1)=0.1;
h=0.25;
n=(tf-t0)/h;
t=t0:h:tf;
for i=1:n
K1=f(t(i),y(i));
K2=f(t(i)+h/2,y(i)+h*K1/2);
K3=f(t(i)+h/2,y(i)+h*K2/2);
K4=f(t(i)+h,y(i)+h*K3);
y(i+1)=y(i)+h*(K1+2*K2+2*K3+K4)/2;
end
plot(t,y,'-*','LineWidth',2)
grid on
title('Runge-Kutta de ordem 4')
Comparação entre os métodos
Funções do matlab para resolver PVIs (site MathWorks)
O ode45 implementa um método Runge-Kutta explícito de ordem intermédia (4,5), conhecido como o par Dormand-Prince
• Runge-Kutta (4,5): Isso significa que ele usa uma combinação de duas fórmulas de Runge-Kutta: uma de quarta ordem e uma de quinta ordem
• Controle de Erro: A diferença entre os resultados dessas duas ordens é usada para estimar o erro de truncamento local
• Passo de Tempo Variável: O ode45 utiliza essa estimativa de erro para ajustar dinamicamente o passo de tempo. Se o erro for muito grande, ele diminui o passo para manter a precisão; se for muito pequeno, ele aumenta o passo para acelerar a computação. Isso o torna muito eficiente para equações diferenciais não rígidas (nonstiff).
Sistemas de EDO: Modelo Presa-Predador
\[\left\{\begin{array}{l} x'(t)=ax-bxy \\ y'(t)=-cy+dxy \end{array}\right.\quad \quad \begin{array}{l} x(t_0)=x_0 \\ y(t_0)=y_0 \end{array} \]
Método de Heun
Aproximação numérica via método de Heun para h igual a 1,0.5,0.25 e 1, e o erro cometido comparado à aproximação via ode45 (Matlab).
Método de Runge-Kutta de 4ª ordem
Aproximação numérica via método de Heun para h igual a 1,0.5,0.25 e 1, e o erro cometido comparado à aproximação via ode45 (Matlab).
RK4 para sistemas
f=@(t,y) [0.25*y(1)-0.01*y(1)*y(2);
-y(2)+0.01*y(1)*y(2)];
%derivadas
t0=0; tf=10;
y(1)=[80; 30];
h=0.25;
n=(tf-t0)/h;
t=t0:h:tf;
for i=1:n
K1=f(t(i),y(:,i));
K2=f(t(i)+h/2,y(:,i)+h*K1/2);
K3=f(t(i)+h/2,y(:,i)+h*K2/2);
K4=f(t(i)+h,y(:,i)+h*K3);
y(:,i+1)=y(:,i)+h*(K1+2*K2+2*K3+K4)/2;
end
plot(t,y,'-*','LineWidth',2)
grid on
title('Runge-Kutta de ordem 4')