Algoritmos Numéricos Fundamentales

Método de Newton-Bisección

f=inline('----');
df=inline('----');
em=eps/2;a=--;b=---;fa=f(a);fb=f(b);x=(a+b)/2;[a x b],fx=f(x);[fa fx fb]

(Si está entre a-x)

b=x;fb=fx;dfx=df(x);if abs(dfx)>em*abs(fx),y=x-fx/dfx;[a y b],else,dfx,end;x=y;fx=f(x);[fa fx fb]

(Si está entre x-b)

a=x;fa=fx;dfx=df(x);if abs(dfx)>em*abs(fx),y=x-fx/dfx;[a y b],else,dfx,end;x=y;fx=f(x);[fa fx fb]

Método de la Secante-Bisección

f=inline('---')
em=eps/2;a=---;b=----;fa=f(a);fb=f(b);x0=a;x1=b;fx0=fa;fx1=fb;m=(fx1-fx0)/(x1-x0);if abs(m)>em*abs(fx1),x=x1-fx1/m;[a x b],else m,end;fx=f(x);[fa fx fb]

(Si está entre x-b)

x0=x1;x1=x;fx0=fx1;fx1=fx;a=x;fa=fx;m=(fx1-fx0)/(x1-x0);if abs(m>em*abs(fx1),x=x1-fx1/m;[a x b],else m,end;fx=f(x);[fa fx fb]

(Si está entre a-x)

x0=x1;x1=x;fx0=fx1;fx1=fx;b=x;fb=fx;m=(fx1-fx0)/(x1-x0);if abs(m>em*abs(fx1),x=x1-fx1/m;[a x b],else m,end;fx=f(x);[fa fx fb]

Deflación de Maehly e itraiz

p=[--------------];  dp=polyder(p);   v=[];  polyval(p,---); polyval(p,---);

% Si encuentras una raíz dando valores

x=---   x=itraiz(p,dp,v,x)

Iteramos

No converge, posible raíz múltiple

Calculamos la raíz del polinomio derivado

d2p=polyder(dp)   u=[];
x=itraiz(dp,d2p,u,x)

%% No converge de manera muy clara

Calculamos la raíz de p''(x)

d3p=polyder(d2p);
x=itraiz(d2p,d3p,u,x)

Comprobamos que es raíz con

polyval(p,x)          polyval(dp,x)            polyval(d2p,x)
v=[v;x;conj(x);x;conj(x);x;conj(x)]    x=i
x=itraiz(p,dp,v,x)      v=[v;x;conj(x)]            polyval(p,x)

Interpolación de Lagrange

x=[----]';
y=[----]';
n=length(x);
v=ones(n,1)
M=[v]
for k=1:n-1;v=v.*x;M=[M v];end
coef=M\y

D es la solución

D=coef(end:-1:1)

Polinomios de Chebyshev en [a,b]

Cálculo de nodos de Chebyshev

n = número de partes en que se divide el intervalo

n=--,a=---,b=---
for k=1:n+1;x(k)=(a+b)/2+((b-a)/2)*cos(2*(n-(k-1))+1)*pi/(2*n+2));end

Interpolación de Hermite

M=[1 a a^2 a^3 ; 1 b b^2 b^3; 0 1 2*a 3*a^2; 0 1 2*b 3*b^2]
b=[-----   ------   ------  ------]'
coef=M\b

Fenómeno de Runge

Ajuste de la función de Runge con polinomio de grado <= 20

En los nodos obtenidos con xj = ..., j = 0, 1, ..., 100. h = 1/20

x=(-5:0.1:5)';
b=inline('1./(1+x.^2)');
b=f(x);
v=ones(length(x),1);
A=[v]
for k=1:20;v=v.*x;A=[A v];end
coef=(A'*A)\(A'*b); % de menor a mayor
k=20
[Q R]=qr(A);
a=R(1:k+1,:)\(Q(:,1:k+1)'*b);
p=polyfit(x,b,20) % de mayor a menor

Integración Numérica Adaptativa

f=inline('---');a=--;b=--;ea=---;er=---;
x=[a (a+b)/2 b]; [i1 e1]=iehardy(x(1),x(2),f);[i2 e2]=iehardy(x(2),x(3),f);I=[i1 i2];E=[e1 e2];itg=sum(I);error=sum(E);[error ea+er*abs(itg)]
[s,k]=max(E)

% Repetir dos líneas, hasta que el error sea menor que ea + er * abs(itg)

% Si tarda mucho

[s,k]=max(E)
ind=1;while ind==1,k1=k+1;xn=(x(k)+x(k1))/2;x=[x(1:k) xn x(k1:end)]; [i1 e2]=iehardy(x(k),xn,f);[i2 e2]=iehardy(xn,x(k1+1),f);
I=[I(1:k-1) i1 i2 I(k1:end)];E=[E(1:k-1) e1 e2 E(k1:end)];itg=sum(I);error=sum(E); if error<ea+er*abs(itg),[itg error],ind=0;else,[s,k]=max(E);end,end

% Para saber el número de iteraciones, usar length(x)