Estadística

Programación con Matlab

Compartir 0 Me sirvió 0 No me sirvió

MÉTODOS DE COMPUTACIÓN ESTADÍSTICA

Chapter 1

Introduction to MATLAB

1.38. The classic quadratic formula says that the two roots of the quadratic equa-

tion

ax2 + bx + c = 0

are

x1; x2 = ¡b § pb2 ¡ 4ac

Use this formula in Matlab to compute both roots for

a = 1; b = ¡100000000; c = 1

Compare your computed results with

roots([a b c])

What happens if you try to compute the roots by hand or with a hand

calculator?

You should ¯nd that the classic formula is good for computing one root, but

not the other. So use it to compute one root accurately and then use the fact

that

x1x2 = c

to compute the other.

X1=(-10^-8+sqrt(10^16-4))/2

X1 =

5.0000e+007

X2=(-10^-8-sqrt(10^16-4))/2

X2 =

-50000000

>> help roots

roots([1 10^8 1])

ans =

1.0e+008 *

-1.0000

-0.0000

Ejercicio sin número sobre error absoluto y error relativo:

Calcular el error absoluto y relativo producido al calcular n factorial con la formula de Stirling n=1,……..,17.

f=factorial(10)

f =

3628800

>> fs=sqrt(2*pi*10)*exp(-10)*(10^10)

fs =

3.598695618741036e+006

>> er=abs(f-fs)/f

er =

0.00829596044394

ea=abs(f-fs)

ea =

3.010438125896407e+004

f=factorial(20)

f =

2.432902008176640e+018

>> fs=sqrt(2*pi*20)*exp(-20)*(20^20)

fs =

2.422786846761134e+018

>> ea=abs(f-fs)

ea =

1.011516141550643e+016

>> er=abs(f-fs)/f

er =

0.00415765262288

Ejercicio sin numero:

% error relativo y absoluto

%en la formula de Stirling

%imprime una tabla mostrando

%el error producido por la

%aproximacion de Stirling

%para n!

close all;clc

disp('')

disp(' Aproximation absolute relative')

disp(' n n! Stirling error error')

disp('----------------------------------------------------------------------------------------------')

e=exp(1);

nfact=1;

for n=1:17

%f=factorial(n);

f=n*f;

fs=sqrt(2*pi*n)*((n/e)^n);

abserror=abs(f-fs);

relerror=abserror/f;

d=sprintf(' %2.0f %16.0f %18.2f %16.2f %5.2e',...

n,f,fs,abserror,relerror);

IMPORTANTE:%18.2 es 18 digitos con 2decimales,por lo tanto quedan 16enteros;la f sigifica el formato que le corresponda y la e significa en formato decimal;como d=sprintf(' %2.0f %16.0f %18.2f %16.2f %5.2e',...

n,f,fs,abserror,relerror) esta entre comillas eso significa que es un texto

disp(d)

end

Lo compruebo o ejecuto:

Aproximation absolute relative

n n! Stirling error error

----------------------------------------------------------------------------------------------

1 1 0.92 0.08 7.79e-002

2 2 1.92 0.08 4.05e-002

3 6 5.84 0.16 2.73e-002

4 24 23.51 0.49 2.06e-002

5 120 118.02 1.98 1.65e-002

6 720 710.08 9.92 1.38e-002

7 5040 4980.40 59.60 1.18e-002

8 40320 39902.40 417.60 1.04e-002

9 362880 359536.87 3343.13 9.21e-003

10 3628800 3598695.62 30104.38 8.30e-003

11 39916800 39615625.05 301174.95 7.55e-003

12 479001600 475687486.47 3314113.53 6.92e-003

13 6227020800 6187239475.19 39781324.81 6.39e-003

14 87178291200 86661001740.60 517289459.40 5.93e-003

15 1307674368000 1300430722199.46 7243645800.54 5.54e-003

16 20922789888000 20814114415223.09 108675472776.91 5.19e-003

17 355687428096000 353948328666099.44 1739099429900.56 4.89e003

14-marzo-06

MATLAB respuestas a la pagina 49 de unidades didácticas 01_intro.pdf

%este programa calcula el

%seno de x usando el

%desarrollo en serie de taylor

%en torno al cero

%sin(x) = x - x^3/3! + x^5/5! - ...

%el criterio de parada del

%bucle "while" consiste en

%que el siguiente termino

%de la serie sea tan pequeño

%que no modifique el resultado

%obtenido hasta el momento

function s = powersin(x)

s = 0;

t = x,

n = 1;

while s + t ~= s;

s = s+t;

t = - x.^2/((n+1)*(n+2)).*t *para hallar ~ es: alt+126

n = n+2;

end

control+c desbloquea las operaciones

lo abrimos y lo usamos:

>> powersin(pi)

t =

3.1416

ans =

2.4791e-016

>> x = 11*(pi)/2;

>> e=abs(powersin(x)-1)

t =

17.2788

e =

2.0000

>> sin(11*pi/2)

ans =

-1

>> sin(21*pi/2)

ans =

>> e = abs(powersin(x)+1)

t =

17.2788

e =

2.1287e-010

*donde e es igual al erros, abs igual al valor absoluto.

>> x = 21*(pi)/2;

>> e = abs(powersin(x)+1)

t =

32.9867

e =

1.9999

>> e = abs(powersin(x)-1)

t =

32.9867

e =

1.3324e-004

>> x = 31*(pi)/2;

>> e = abs(powersin(x)-1)

t =

48.6947

e =

5.8230e+003

modificamos el programa inicial:

function [s,k] = powersin(x)

s = 0;

t = x,

n = 1;

k = 0;

while s + t ~= s;

k = k +1;

s = s+t;

t = - x.^2/((n+1)*(n+2)).*t;

n = n+2;

end

Operamos:

>> [s,k]= powersin(x)

t =

48.6947

s =

-5.8220e+003

k =

>> x = (pi)/2;

>> [s,k]= powersin(x)

t =

1.5708

s =

1.0000

k =

*cambiamos de nuevo el programa :

function [s,k,tmax] = powersin(x)

s = 0;

t = x,

n = 1;

k = 0;

tmax=x;

while s + t ~= s;

k = k +1;

s = s+t;

t = - x.^2/((n+1)*(n+2)).*t;

if t > tmax

tmax = t;

end

n = n+2;

end

operamos:

>> x = (pi)/2;

>> [s,k,tmax]= powersin(x)

t =

1.5708

s =

1.0000

k =

tmax =

1.5708

>> x = 11*(pi)/2;

>> [s,k,tmax]= powersin(x)

t =

17.2788

s =

-1.0000

k =

tmax =

3.0665e+006

Modificando de nuevo:

function [s,k,tmax,t] = powersin(x)

s = 0;

t = x,

n = 1;

k = 0;

tmax=x;

while s + t ~= s;

k = k +1;

s = s+t;

t = - x.^2/((n+1)*(n+2)).*t;

if t > tmax

tmax = t;

end

n = n+2;

end

Operamos

>> [s,k,tmax,t]= powersin(x)

t =

17.2788

s =

-1.0000

k =

tmax =

3.0665e+006

t =

-2.6232e-017

28-marzo-06

Tenemos el programa lutx:

function [L,U,p] = lutx(A)

%LUTX Triangular factorization, textbook version

% [L,U,p] = lutx(A) produces a unit lower triangular matrix L,

% an upper triangular matrix U, and a permutation vector p,

% so that L*U = A(p,:)

[n,n] = size(A);

p = (1:n)';

for k = 1:n-1

% Find index of largest element below diagonal in k-th column

[r,m] = max(abs(A(k:n,k)));

m = m+k-1;

% Skip elimination if column is zero

if (A(m,k) ~= 0)

% Swap pivot row

if (m ~= k)

A([k m],:) = A([m k],:);

p([k m]) = p([m k]);

end

% Compute multipliers

i = k+1:n;

A(i,k) = A(i,k)/A(k,k);

% Update the remainder of the matrix

j = k+1:n;

A(i,j) = A(i,j) - A(i,k)*A(k,j);

end

y aplicamos en matlab:

>>A= [0 1;1 0]

A =

0 1

1 0

>> [L,U,p]= lutx(A)

function [L,U,p,sig] = lutx(A)

%LUTX Triangular factorization, textbook version

% [L,U,p] = lutx(A) produces a unit lower triangular matrix L,

% an upper triangular matrix U, and a permutation vector p,

% so that L*U = A(p,:)

[n,n] = size(A);

p = (1:n)';

sig=1;

for k = 1:n-1

% Find index of largest element below diagonal in k-th column

[r,m] = max(abs(A(k:n,k)));

m = m+k-1;

% Skip elimination if column is zero

if (A(m,k) ~= 0)

% Swap pivot row

if (m ~= k)

A([k m],:) = A([m k],:);

p([k m]) = p([m k]);

sig = -1*sig;

end

% Compute multipliers

i = k+1:n;

A(i,k) = A(i,k)/A(k,k);

% Update the remainder of the matrix

j = k+1:n;

A(i,j) = A(i,j) - A(i,k)*A(k,j);

end

% Separate result

L = tril(A,-1) + eye(n,n);

U = triu(A);

Aplicamos::

>> A

A =

0 1

1 0

>> [L,U,p] = lutx(A)

L =

1 0

0 1

U =

1 0

0 1

p =

>> [L,U,p,sig] = lusigno(A)

> [L,U,p,sig] = lusigno(A)

L =

1 0

0 1

U =

1 0

0 1

p =

sig =

-1

>> A = [1 2 3;2 1 4;-1 2 1]

A =

1 2 3

2 1 4

-1 2 1

>> [L,U,p,sig] = lusigno(A)

L =

1.0000 0 0

-0.5000 1.0000 0

0.5000 0.6000 1.0000

U =

2.0000 1.0000 4.0000

0 2.5000 3.0000

0 0 -0.8000

p =

sig =

Programa:

%Ejercicio 2.7 del libro Moler

%este programa calcula el determinante

%de la matriz A utilizando la factorización

%LU como de A ,PA=LU

%det(A) = +- U(1,1)*U(2,2)*...*U(n,n)

%donde U(i,i) son los elementos de la diagonal

%de la matriz U, y el signo +- es

%+ si el numero de intercambios de filas

%que hemos realizado es par y menos en caso

%contrario. El signo se almacena en la

%variable sig

function d = mideterminante(A)

[L,U,p,sig] = lusigno(A)

Aplico:

>> A

A =

1 2 3

2 1 4

-1 2 1

>> A = [0 1;1 0]

A =

0 1

1 0

>> diag(A)

ans =

>> X

>> X=diag(A)

X =

>> B=diag(X)

B =

0 0

>> X= 1:5

X =

1 2 3 4 5

B =

1 0 0 0 0

0 2 0 0 0

0 0 3 0 0

0 0 0 4 0

0 0 0 0 5

>> mideterminante(A)

ans =

-1

> B=eye(50)

B =

Columns 1 through 12

1 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0 0 0

Columns 13 through 24

0 0 0 0 0 0 0 0 0 0 0 0

1 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0 0 0

Columns 25 through 36

0 0 0 0 0 0 0 0 0 0 0 0

1 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0 0 0

Columns 37 through 48

0 0 0 0 0 0 0 0 0 0 0 0

1 0 0 0 0 0 0 0 0 0 0 0

0 1 0 0 0 0 0 0 0 0 0 0

0 0 1 0 0 0 0 0 0 0 0 0

0 0 0 1 0 0 0 0 0 0 0 0

0 0 0 0 1 0 0 0 0 0 0 0

0 0 0 0 0 1 0 0 0 0 0 0

0 0 0 0 0 0 1 0 0 0 0 0

0 0 0 0 0 0 0 1 0 0 0 0

0 0 0 0 0 0 0 0 1 0 0 0

0 0 0 0 0 0 0 0 0 1 0 0

0 0 0 0 0 0 0 0 0 0 1 0

0 0 0 0 0 0 0 0 0 0 0 1

0 0 0 0 0 0 0 0 0 0 0 0

Columns 49 through 50

0 0

1 0

0 1

Pasamos al ejemplo 2.11 del libro de moler capitulo 002 pagina 37.

2.7 lutx, bslashtx, lugui

We have three functions implementing the algorithms discussed in this chapter.

The ¯rst function, lutx, is a readable version of the built-in Matlab function lu.

There is one outer for loop on k that counts the elimination steps. The inner loops

on i and j are implemented with vector and matrix operations, so that the overall

function is reasonably e±cient.

function [L,U,p] = lutx(A)

%LU Triangular factorization

% [L,U,p] = lutx(A) produces a unit lower triangular

% matrix L, an upper triangular matrix U, and a

% permutation vector p, so that L*U = A(p,:).

[n,n] = size(A);

p = (1:n)'

for k = 1:n-1

% Find largest element below diagonal in k-th column

[r,m] = max(abs(A(k:n,k)));

m = m+k-1;

% Skip elimination if column is zero

if (A(m,k) ~= 0)

% Swap pivot row

if (m ~= k)

A([k m],:) = A([m k],:);

p([k m]) = p([m k]);

end

% Compute multipliers

i = k+1:n;

A(i,k) = A(i,k)/A(k,k);

% Update the remainder of the matrix

j = k+1:n;

A(i,j) = A(i,j) - A(i,k)*A(k,j);

end

% Separate result

L = tril(A,-1) + eye(n,n);

U = triu(A);

Y la versión simplificada de este operador la llama: Bslashtx(A,b) [backslash]:

function x = bslashtx(A,b)

% BSLASHTX Solve linear system (backslash)

% x = bslashtx(A,b) solves A*x = b

[n,n] = size(A);

if isequal(triu(A,1),zeros(n,n))

% Lower triangular

x = forward(A,b);

return

elseif isequal(tril(A,-1),zeros(n,n))

% Upper triangular

x = backsubs(A,b);

return

elseif isequal(A,A')

[R,fail] = chol(A);

if ~fail

% Positive definite

y = forward(R',b);

x = backsubs(R,y);

return

end

Realizamos otro programa para la inversa (ejercicio 2.11 de moler)

%Ejercicio 2.11 del libro de Moler

%este programa nos permite calcular

%la matriz inversa de una matriz

%usando la factorización LU y

%las versiones sencillas de los

%programas LU y \ que nos ofrece

%el libro de Moler

function X = miinversa(A)

[n,n] = size (A);

E = eye(n);

for i=1:n

%x(:,i)=bslashtx(A,E(:,i));

X(:,i)= A\E(:,i);

end

aplicamos en matlab:

A =

0 1

1 0

>> miinversa(A)

ans =

0 1

1 0

>> B=rand(50);

>> X=miinversa(B);

>> norm(B*X - eye(50))

ans =

2.3934e-014

>> X2=inv(B);

>> norm(B*X2 - eye(50))

ans =

1.2760e-013

>> cond(B)

ans =

517.6219

4-abril-06

CAPITULO 3 PDF

A =

1 1 1

1 2 4

1 3 9

>> B= [2;1;3]

B =

>> P = A\B

P =

6.0000

-5.5000

1.5000

>> X= [1;2;3]

X =

>> Y = [2; 1; 3]

Y =

>>plot(X,Y,'r*')

Y obtenemos un gráfico:

>> plot(X,Y,'r*',u,v)

u=0:0.01:4;

>> v= 6 -5.5*u+315*u.^2;

(es el polinomio obtenido antes evaluado en los puntos dados)

>> plot(X,Y,'r*',u,v)

Y obtenemos:

*Tenemos los siguientes programas:

function v = polyinterp(x,y,u)

%POLYINTERP Polynomial interpolation.

% v = POLYINTERP(x,y,u) computes v(j) = P(u(j)) where P is the

% polynomial of degree d = length(x)-1 with P(x(i)) = y(i).

% Use Lagrangian representation.

% Evaluate at all elements of u simultaneously.

n = length(x);

v = zeros(size(u));

for k = 1:n

w = ones(size(u));

for j = [1:k-1 k+1:n]

w = (u-x(j))./(x(k)-x(j)).*w;

end

v = v + w*y(k);

end

número de puntos

function v = piecelin(x,y,u)

%PIECELIN Piecewise linear interpolation.

% v = piecelin(x,y,u) finds the piecewise linear L(x)

% with L(x(j)) = y(j) and returns v(k) = L(u(k)).

% First divided difference

delta = diff(y)./diff(x);

% Find subinterval indices k so that x(k) <= u < x(k+1)

n = length(x);

k = ones(size(u));

for j = 2:n-1

k(x(j) <= u) = j;

end

% Evaluate interpolant

s = u - x(k);

v = y(k) + s.*delta(k);

>> X= 0:10

X =

0 1 2 3 4 5 6 7 8 9 10

>> Y=rand(1,11);

>> u= -0.1:0.01:10.1;

>> Y = 10*Y;

>> v= polyinterp(X,Y,u);

>> plot(X,Y,'r*',u,v)

Y obtenemos:

>> x=[1 4 5 4 -1 2 0]

x =

1 4 5 4 -1 2 0

>> diff(x)

ans =

3 1 -1 -5 3 -2

>> x

x =

0 1 2 3 4 5 6 7 8

>> y = 10*rand(1,9);

>> y

y =

Columns 1 through 7

4.3866 4.9831 2.1396 6.4349 3.2004 9.6010 7.2663

Columns 8 through 9

4.1195 7.4457

>> u = -1:0.01:9;

>> v = piecelin(x,y,u);

>> plot(x,y,'r*',u,v)

06/04/06

Del capitulo 3pdf la segunda parte.

'Programación con Matlab'

Este es un buen método (el pichip) para interpolar, pero los hay mejores(el splines)

*Los ptos di son “knots”

K=1 d0+4d1+d3 = 31+3 2 eliminada

K=2 d1+4d2+d4=3 2+3 3

K=3 d2+4d3+d5=3 3+3 4

K=4 d1+4d4+d6=3 4+3 5

K=5 d4+4d5+ d7=3 5+3 6 eliminada

.-Si esto no funciona usamos spline, de la cual una versión simplificada es la splinetx.m , que es:

function v = splinetx(x,y,u)

%SPLINETX Textbook spline function.

% v = splinetx(x,y,u) finds the piecewise cubic interpolatory

% spline S(x), with S(x(j)) = y(j), and returns v(k) = S(u(k)).

% See SPLINE, PCHIPTX.

% First derivatives

h = diff(x);

delta = diff(y)./h;

d = splineslopes(h,delta);

% Piecewise polynomial coefficients

n = length(x);

c = (3*delta - 2*d(1:n-1) - d(2:n))./h;

b = (d(1:n-1) - 2*delta + d(2:n))./h.^2;

% Find subinterval indices k so that x(k) <= u < x(k+1)

k = ones(size(u));

for j = 2:n-1

k(x(j) <= u) = j;

end

% Evaluate spline

s = u - x(k);

v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k)));

Ensaya las primeras derivadas

* Implementa alguno de los métodos para hallar las di. (not-a kan ; loked)

De momento esto lo olvidamos

P(x) = s, es el output.

*para ordenar los puntos de manera creciente, en matlab usamos el comando “sort”.

Por ejemplo si tenemos un vector X usamos un comando Y=sort(x), nos aparece un vector ordenado de menor a mayor.

miramos del capitulo 3:

'Programación con Matlab'

Donde gráficamente reconocemos si una función es o no derivable por las propiedades estudiadas en cálculo.

Seguimos usando el programa splinetx.m y aplicamos en matlab.

*el punto antes de una operación significa, componente a componente.

También usamos el programa splineslopes:

function d = splineslopes(h,delta)

% SPLINESLOPES Slopes for cubic spline interpolation.

% splineslopes(h,delta) computes d(k) = S'(x(k)).

% Uses not-a-knot end conditions.

% Diagonals of tridiagonal system

n = length(h)+1;

a = zeros(size(h)); b = a; c = a; r = a;

a(1:n-2) = h(2:n-1);

a(n-1) = h(n-2)+h(n-1);

b(1) = h(2);

b(2:n-1) = 2*(h(2:n-1)+h(1:n-2));

b(n) = h(n-2);

c(1) = h(1)+h(2);

c(2:n-1) = h(1:n-2);

% Right-hand side

r(1) = ((h(1)+2*c(1))*h(2)*delta(1)+h(1)^2*delta(2))/c(1);

r(2:n-1) = 3*(h(2:n-1).*delta(1:n-2)+h(1:n-2).*delta(2:n-1));

r(n) = (h(n-1)^2*delta(n-2)+(2*a(n-1)+h(n-1))*h(n-2)*delta(n-1))/a(n-1);

% Solve tridiagonal linear system

d = tridisolve(a,b,c,r);

Usa una matriz tri-diagonal (con tres diagonales) y se soluciona gracias a tridisolve

Tridiagonal:

'Programación con Matlab'

Trisolve(a,b,c,r)

Ax = r

x = A \ r

las condiciones descritas en nuestro ejemplo “not-a-knot”

'Programación con Matlab'

Para la primera ecuación: h2d1+h2+a1d2 = … = (hn-1+hn-2)dn-1+hn-2*dn= …

Del capitulo 3 del libro de Mooler, página 22 :

3.11. Modify splinetx and pchiptx so that, if called with two output arguments,

they produce both the value of the interpolant and its ¯rst derivative. That

is,

[v,vprime] = pchiptx(x,y,u)

and

[v,vprime] = splinetx(x,y,u)

compute P(u) and P0(u).

[v,v1,v2,v3] = splintex4(x,y,n)

Creamos un nuevo programa a partir del splintex:

%Este programa resuelve el problema 3.11

% del libro de mooler(modificado)

function [v,v1,v2,v3] = splinetx3(x,y,u)

%SPLINETX Textbook spline function.

% v = splinetx(x,y,u) finds the piecewise cubic interpolatory

% spline S(x), with S(x(j)) = y(j), and returns v(k) = S(u(k)).

% and the first, second and third derivatives

%v1 , v2, v3.

% See SPLINE, PCHIPTX.

% First derivatives

h = diff(x);

delta = diff(y)./h;

d = splineslopes(h,delta);

% Piecewise polynomial coefficients

n = length(x);

c = (3*delta - 2*d(1:n-1) - d(2:n))./h;

b = (d(1:n-1) - 2*delta + d(2:n))./h.^2;

% Find subinterval indices k so that x(k) <= u < x(k+1)

k = ones(size(u));

for j = 2:n-1

k(x(j) <= u) = j;

end

% Evaluate spline and its three first derivatives.

s = u - x(k);

v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k)));

v1= d(k) +s.*(2*c(k)) + 3*s.*(b(k));

v2= 2+c(k)+6*s.*b(k);

v3= 6*b(k);

%-----------------------------------------------------------------------

En command window:

>> x = [1 1.5 3 4.5 7 8.1 9.2 10]

x =

Columns 1 through 7

1.0000 1.5000 3.0000 4.5000 7.0000 8.1000 9.2000

Column 8

10.0000

>> y = [2 -3 4 6 2.3 4 -1 3 6 8]

y =

Columns 1 through 7

2.0000 -3.0000 4.0000 6.0000 2.3000 4.0000 -1.0000

Columns 8 through 10

3.0000 6.0000 8.0000

>> format short

>> x

x =

Columns 1 through 7

1.0000 1.5000 3.0000 4.5000 7.0000 8.1000 9.2000

Column 8

10.0000

>> u = 0.5:0.01:10.5;

>> [v,v1,v2,v3]=splinetx3(x,y,u)

>> [v,v1,v2,v3]=splinetx3(x,y,u);

>> plot(x,y,'r*',u,v2,'g',u,v3,'c')

>> plot(x,y,'r*',u,v,'r*',u,v1,'b',u,v2,'g',u,v3,'c')

Teniamos que tener el programa triidsolve:

function x = tridisolve(a,b,c,d)

% TRIDISOLVE Solve tridiagonal system of equations.

% x = TRIDISOLVE(a,b,c,d) solves the system of linear equations

% b(1)*x(1) + c(1)*x(2) = d(1),

% a(j-1)*x(j-1) + b(j)*x(j) + c(j)*x(j+1) = d(j), j = 2:n-1,

% a(n-1)*x(n-1) + b(n)*x(n) = d(n).

% The algorithm does not use pivoting, so the results might

% be inaccurate if abs(b) is much smaller than abs(a)+abs(c).

% More robust, but slower, alternatives with pivoting are:

% x = T\d where T = diag(a,-1) + diag(b,0) + diag(c,1)

% x = S\d where S = spdiags([[a; 0] b [0; c]],[-1 0 1],n,n)

x = d;

n = length(x);

for j = 1:n-1

mu = a(j)/b(j);

b(j+1) = b(j+1) - mu*c(j);

x(j+1) = x(j+1) - mu*x(j);

end

x(n) = x(n)/b(n);

for j = n-1:-1:1

x(j) = (x(j)-c(j)*x(j+1))/b(j);

end

cambiamos le valor de y:

>> y = [0 2 4 -1 -3 0 4 8 9 9];

>> u=0.9:0.01:10.1;

>> [v,v1,v2,v3] = splinetx3(x,y,u);

>> plot(x,y,'r*',u,v,'r*',u,v1,'b',u,v2,'g',u,v3,'c')

Mirar problem3.11

Mirar problem 3.9

18-ABRIL-06

<<Chapter 4>>

Teorema de los ceros de Bolzano.

*en matlab el poner ó es lo mismo que poner =.

En los apuntes:

Given a,b,f(a),f(b),e

while |b-a|>e|b|

x=a+b/2

if f(x) = 0 then

end

if f(x)f(a)<0 then

b=x

else

a=x

end

Hacemos un programa para hala ceros con el método de la bisección:

%este programa calcula los ceros

%de una función f(x) por el

%método de la bisección

%Inputs a,b,e

%a,b extremos del intervalo

%e error relativo admisible

% la función f(x) se escribe

%en un fichero auxiliar de nombre

%fun.m

function c = bisec(a,b,e)

fa = feval(`fun',a);

fb = feval(`fun',b);

while abs(b-a) > e*abs(b)

c = (a+b)/2;

fc = feval(`fun',c);

if fc == 0

break

end

if fc*fa < 0

b = c;

else

a = c;

end

y otro programa:

% este fichero sirve para almacenar la función

%f(x)

function y = fun(x)

y = sin(x.^2)-x.^3;

En matlab:

>> c = feval('fun',0)

c =

>> feval('fun',3)

ans =

-26.5879

>>>> sin(3^2)-3^3

ans =

-26.5879

>>>> feval('fun',8)

ans =

-511.0800

>> feval('fun',0.1)

ans =

0.0090

>> feval('fun',0)

ans =

>> c= bisec(0.1,3,10^-6)

c =

0.8960

>> format long

>> c= bisec(0.1,3,10^-6)

c =

0.89599368572235

Modificamos el primer programa:

%este programa calcula los ceros

%de una función f(x) por el

%método de la bisección

%Inputs a,b,e

%a,b extremos del intervalo

%e error relativo admisible

% la función f(x) se escribe

%en un fichero auxiliar de nombre

%fun.m

function [c,k] = bisec(a,b,e)

k=0;

fa = feval('fun',a);

fb = feval('fun',b);

while abs(b-a) > e*abs(b)

k = k+1;

c = (a+b)/2;

fc = feval('fun',c);

if fc == 0

break

end

if fc*fa < 0

b = c;

else

a = c;

end

modificamos también el segundo:

% este fichero sirve para almacenar la función

%f(x)

function y = fun(x)

y = x.^3 - 2.000001*x.^2+1.000002*x - 0.000001;

>> feval('fun',-10)

ans =

-1.210000121000000e+003

>> feval('fun',10)

ans =

8.099999190000000e+002

>> c = bisec(-10,10,10^-6)

c =

1.000000224848918e-006

>> u = -2:0.01:2;

>> v = feval('fun',u);

>> plot(u,v)

>> z= 0*u;

>> plot(u,v,u,z,'k')

>> c = bisec(-1,0.5,10^-12)

c =

1.000000000000350e-006

>> c = bisec(-1,0.5,10^-13)

c =

9.999999999999430e-007

>> c = bisec(-1,0.5,10^-14)

c =

9.999999999999989e-007

>> c = bisec(-1,0.5,10^-15)

c =

9.999999999999995e-007

>> c = bisec(-1,0.5,10^-16)

c =

1.000000000000000e-006

>> u = 0.5:0.001:1.5;

>> z= 0*u;

>> v = feval('fun',u);

>> plot(u,v,u,z,'k')

>> feval('fun',1)

ans =

-8.226659269441623e-017

Otro método es: EL MÉTODO DE NEWTON

Creamos el programa newton:

%Este programa obtiene los

%ceros de la función f(x)

% por el método de Newton

%Inputs: x0 (punto de partida),e (error relativo)

%Necesitamos dos ficheros auxiliares:

% - fun.m (donde guardamos la función f(x))

% - derfun.m (donde almacenamos la derivada

%f'(x) de la función f(x))

function [c,k] = newton(a,e)

k=0;

b=a;

fa=feval('fun',a);

dfa=feval('derfun',a);

c=a-fa/dfa;

while abs(c-b) > e*abs(b)

a=c;

b=a;

k = k+1;

fa = feval ('fun',a);

dfa=feval('derfun',a);

c= a-fa/dfa;

end

creamos el programa derfun:

%este fichero contiene la derivada de la función

%f(x) que esta guardada en el fichero fun.m

function dy = derfun(x)

%la otra función que usabamos era: dy = 2*cos (x.^2).*x-3*x.^2;

dy = 3*x^2-2*2.000001*x+1.000002;

20-abril-06

+página 10 del capítulo 5

Fix point method

(puntos fijos)

'Programación con Matlab'

El algoritmo de punto fijo dice:

Dados x y e mientras que |g(x) -x | > |x|:

Given x

While |g(x) -x | > |x|

X g(x)

End

***** Creamos un programa llamado fijo.m: *******

%este programa calcula los puntos

%fijos de la función g(x)

% g(x*) = x*

% la función g(x) esta almacenada en el fichero fun.m

%inputs: dato inicio x0, error relativo (epsilon)

% outputs: punto fijo x*, numero de iteraciones k

function [x,k] = fijo(x0,e)

k = 0;

x= x0;

gx = feval('fun',x);

% en feval('fun',x); evalua g(x)

while abs( gx - x) > e*abs(x)

x = gx;

gx=feval('fun',x);

k = k+1;

end

y el programa:

function y = fun(x)

y = cos(x);

y ejecuto en matlab:

>> [x,k]=fijo (0,10^-6)

x =

0.7391

k =

Y con el programa derfun Nuevo:

%este fichero contiene la derivada de la función

%f(x) que esta guardada en el fichero fun.m

function dy = derfun(x)

%la otra función que usabamos era: dy = 2*cos (x.^2).*x-3*x.^2;

%dy = 3*x^2-2*2.000001*x+1.000002;

dy = -sin(x)-1;

en matlab:

>> [x,k] = newton(0.1,10^-6)

x =

0.7391

k =

'Programación con Matlab'

26-abril-06

+apuntes libro moler pagina zero página 20 en adelante

Con los programas:

derfun:

function dy = dfun(x)

% y = cos(x);

dy = 3*x^2-2 -5;

fun:

function y = fun(x)

%y = cos(x)-x;

y = x.^3-2.*x-5;

En matlab:

>> u = 0:0.01:3;

>> v = fun(u);

>> z=0*u;

>> plot(u,v,u,z,'k');

>> [c,k] = bisec(0,3,eps)

c =

2.0946

k =

>> [c,k] = newton(0,eps)

c =

2.0946

k =

Donde bisec:

%Este programa calcula los ceros de una funcion f(x)

%por el metodo de la biseccion

%e error relativo admisible

%la funcion f(x) se escribe en un fichero auxiliar

%con el nombre fun.m

function [x,p] = bisec(a,b,e)

fa = feval('fun',a); %feval evalua la funcion en el punto dado

fb = feval('fun',b);

p = 0;

while abs(b-a) > e*abs(b)

p = p+1;

x = (a+b)/2;

fx = feval('fun',x);

if fx == 0

break

end

if fx*fa < 0

b = x;

else a = x;

end

Donde newton:

%Este programa obtiene los ceros de la funcion f(x)

%por el metodo de Newton

%x0 = punto de partida, e = error relativo

%necesitamos dos ficheros auxiliares

%fun.m(donde guardamos la funcion f(x)

%derfun.m(donde almacenamos la derivada

%f'(x) de la funcion f(x)

function [c,k] = newton(x,e)

k = 0;

b = x;

fx = feval('fun',x);

dfx = feval('derfun',x);

c = x - fx/dfx;

while abs(c-b) > e*abs(b);

x = c;

b = x;

k = k+1;

fx = feval('fun',x);

dfx = feval('derfun',x);

c = x-fx/dfx;

end

la nueva function derfun:

function dy = derfun(x)

% y = cos(x);

dy = 1./(1+x.^2);

>> [c,k] = newton(0,eps)

c =

1.7321

k =

Ahora la funcion fun:

function y = fun(x)

y = sign(x-2).*sqrt(abs(x-2));

>> v = fun(u);

>> u = 0:0.01:4;

>> v = fun(u);

>> z=0*u;

>> plot(u,v,u,z,'k');

Ahora derfun:

function dy = derfun(x)

% y = cos(x);

%dy = 1./(1+x.^2);

dy = 1./(2*sqrt(abs(x-2)));

>> [c,k] = bisec(0,4,eps)

c =

k =

>> v = fun(u);

>> u = 0:0.01:4;

>> v = fun(u);

>> z=0*u;

>> plot(u,v,u,z,'k');

>> [c,k] = bisec(0,4,eps)

c =

k =

[c,k] = newton (0,eps)

'Programación con Matlab'

>> u = -20:0.01:20;

>> v = besselj(0,u);

>> z=0*u;

>> plot(u,v,u,z,'k');

Warning: Imaginary parts of complex X and/or Y arguments ignored.

>> u = 0:0.01:20;

>> v = besselj(0,u);

>> z=0*u;

>> plot(u,v,u,z,'k');

>> u = 0:0.01:40;

>> v = besselj(0,u);

>> z=0*u;

>> plot(u,v,u,z,'k');

>> u = -69:0.01:69;

>> v = besselj(0,u);

>> z=0*u;

>> plot(u,v,u,z,'k');

Cambio fun:

function y = fun(x)

%y = cos(x)-x;

%y = x.^3-2.*x-5;

%y = atan(x)-pi/3;

%y = sign(x-2).*sqrt(abs(x-2));

y = besselj(0,x);

y creo el problema 4 el programa:

% este pequeño programa (script)

%resuelve el problema 4.10

% del libro de mooler, buscamos

%los 10 primeros ceros de la función

% de Bessel de primera especie

%de orden 0, J_0(x),

%Para ello usamos el metodo de la bisección

%aplicado a los intervalos [1,4],

% [4,7],...,[29,32]

%Los intervalos los hemos obtenido

%por inspección representando graficamente

%la función.

x = [1 4 7 10 13 17 20 24 27 29 32];

c = zeros(10,1);

for n = 1:10

c(n) = bisec(x(n),x(n+1),eps)

end

*que no lleva function porque no tiene inputs

>> problem4_10

c =

2.4048

c =

2.4048

5.5201

c =

2.4048

5.5201

8.6537

c =

2.4048

5.5201

8.6537

11.7915

c =

2.4048

5.5201

8.6537

11.7915

14.9309

c =

2.4048

5.5201

8.6537

11.7915

14.9309

18.0711

c =

2.4048

5.5201

8.6537

11.7915

14.9309

18.0711

21.2116

c =

2.4048

5.5201

8.6537

11.7915

14.9309

18.0711

21.2116

24.3525

c =

2.4048

5.5201

8.6537

11.7915

14.9309

18.0711

21.2116

24.3525

27.4935

c =

2.4048

5.5201

8.6537

11.7915

14.9309

18.0711

21.2116

24.3525

27.4935

30.6346

>> besselj(0,c(1))

ans =

09/05/06

Capitulo 5 Last squares (minimos cuadrados)

Recordamos:

Reales(n) Complejos(n)

( , ) < , >

Linealidad:

< x, y+z > = <x,y> + <x,z>

<x,By> = B <x,y>

<x,x> mayor o igual que cero, <x,y> + <x,z>

<x,y> = 1/ <x,y>

(x,y) = xT * y ; <x,y> = media(xT) * y

[ e1,e2,…,en] ortonormal

<ei,ej> = dij = 1 si i = j; 0 si i distinto de j

En matlab:

>> Y = [-1,3,0];

>> X=[1,2,3];

>> X*Y'

ans =

>> X=[i,-1,0]

X =

0 + 1.0000i -1.0000

X =

0 + 1.0000i -1.0000 0

>> Y = [-2; sqrt(2) + i; 1]

Y =

-2.0000

1.4142 + 1.0000i

1.0000

y =

-2.0000

1.4142 + 1.0000i

1.0000

>> Y'*Y

ans =

>> X= [i;-1;0]

X =

0 + 1.0000i

-1.0000

>> X'*X

ans =

>> X'*Y

ans =

-1.4142 + 1.0000i

>> norm(X)

ans =

1.4142 (es raís de2)

>> norm(Y)

ans =

(es raíz de8)

*tenemos un vector fila: (v1,….,vn) no ortonormal

aplicamos Grand Smith (q1,…,qn) ortonormal

<qi,qj> = dij

Peo no nos sirve Grand Smith, asi que usamos otros métodos. En matlab tenemos la factorización de QR.

Es decir: tenemos = |v1|, … |vn| qr Q 0 |q1|,…|qn|

En matlab:

>> A= rand(4)

A =

0.9501 0.8913 0.8214 0.9218

0.2311 0.7621 0.4447 0.7382

0.6068 0.4565 0.6154 0.1763

0.4860 0.0185 0.7919 0.4057

>> rank(A)

ans =

>>> v1 = A(:,1)

v1 =

0.9501

0.2311

0.6068

0.4860

>> v2 = A(:,2)

v2 =

0.8913

0.7621

0.4565

0.0185

>> v3 = A(:,3)

v3 =

0.8214

0.4447

0.6154

0.7919

>> v4 = A(:,4)

v4 =

0.9218

0.7382

0.1763

0.4057

>> v2'*v3

ans =

1.3666

>> v1'*v4

ans =

1.3506

>> [Q,R]= qr(A)

Q =

-0.7606 -0.1354 0.4523 0.4457

-0.1850 -0.8151 -0.5489 -0.0062

-0.4858 0.0754 0.0617 -0.8686

-0.3890 0.5582 -0.7002 0.2163

R =

-1.2492 -1.0478 -1.3141 -1.0811

0 -0.6971 0.0148 -0.4867

0 0 -0.3892 -0.2615

0 0 0 0.3409

q3 =

0.4523

-0.5489

0.0617

-0.7002

>> q4=Q(:,4)

q4 =

0.4457

-0.0062

-0.8686

0.2163

>> q1 =

-0.7606

-0.1850

-0.4858

-0.3890

>> q2=Q(:,2)

q2 =

-0.1354

-0.8151

0.0754

0.5582

>> q2'*q3

ans =

5.5511e-017

>> q1'*q4

ans =

1.3878e-017

>> q1'*q1

ans =

1.0000

>> q4'*q4

ans =

1.0000

*perpendicular X e Y : X_|_ Y

En los apuntes:

Classical Gram-Schmidt algorithm:

for j = 1:n

vj = aj

for k = 1:j-1

rjk = aj*qk

vj= vj - rjk qk

end

rjj = ||vj||

qj = vj/rjj

end

Exercise: Study the numerical stability of the

classical GS algorithm.

Exercise: Modify the classical GS algorithm

to correct its numerical unstability.

Gram-Schmidt

orthogonallization method:

Initial basis B = {ai }

ui ð

, ri

= vi ,vi

vi ð ai ð ai ,uk uk

kð1

ið1

Classical Gram-Schmidt algorithm:

for j = 1:n

vj = aj

for k = 1:j-1

rjk = aj*qk

vj= vj - rjk qk

end

rjj = ||vj||

qj = vj/rjj

end

Exercise: Study the numerical stability of the

classical GS algorithm.

Exercise: Modify the classical GS algorithm

to correct its numerical unstability.

QR factorization:

A = QR

Q = [u1 u2 … un] orthogonal matrix

R = [Rij] upper triangular matrix

Tengo los siguientes programas:

Demo_ortog:

% This script explores the lost of orthogonality

% due to the intrinsic numerical unstability

% of the classical and

% modified Gram-Schmidt algorithms.

clc

disp(' Losing orthogonality on Gram-Schmidt algorithms ')

disp(' delta Error on QR Error on CGS Error on MGS ')

disp('------------------------------------------------------------------------')

for delta = logspace(-4,-16,7)

A = [1+delta 2;1 2];

[Q1,R1]= qr(A);

[Q2,R2]= cgs(A);

[Q3,R3]= mgs(A);

error1 = norm(Q1'*Q1 - eye(2));

error2 = norm(Q2'*Q2 - eye(2));

error3 = norm(Q3'*Q3 - eye(2));

disp(sprintf(' %5.0e %20.15f %20.15f %20.15f ', delta,error1,error2,error3))

end

cgs:

% Gram-Schmidt orthogonalization (unstable)

function [Q,R] = cgs(A)

n = length(A);

% starting variables

R = zeros(n,n);

Q = zeros(n,n);

V = zeros(n,n);

% Normalizing the first column vector A

R(1,1)=norm(A(:,1));

Q(:,1)= A(:,1)/R(1,1);

% Classical Gram-Schmidt

for j = 2:n

V(:,j) = A(:,j);

for i = 1:j-1

R(i,j) = Q(:,i)'*A(:,j);

V(:,j)= V(:,j)-R(i,j)*Q(:,i);

end

R(j,j)=norm(V(:,j));

Q(:,j)=V(:,j)/R(j,j);

end

Mgs:

% Modified Gram-Schmidt algorithm (stable)

function [Q,R] = mgs(A)

[m,n] = size(A);

% Initializating variables

R = zeros(m,n);

S = zeros(m,m);

Q = zeros(m,m);

V = [A floor(10*rand(m,m-n))];

for i = 1:n

R(i,i)=norm(V(:,i));

Q(:,i)=V(:,i)/R(i,i);

for j = i+1:n

R(i,j) = Q(:,i)'*V(:,j);

V(:,j)= V(:,j)-R(i,j)*Q(:,i);

end

for i = n+1:m

S(i,i)=norm(V(:,i));

Q(:,i)=V(:,i)/S(i,i);

for j = i+1:n

S(i,j) = Q(:,i)'*V(:,j);

V(:,j)= V(:,j)-S(i,j)*Q(:,i);

end

En matlab:

Demo_ortog

Losing orthogonality on Gram-Schmidt algorithms

delta Error on QR Error on CGS Error on MGS

------------------------------------------------------------------------

1e-004 0.000000000000000 0.000000000002109 0.000000000002109

1e-006 0.000000000000001 0.000000000399139 0.000000000399139

1e-008 0.000000000000000 0.000000028306691 0.000000028306691

1e-010 0.000000000000000 0.000003330622496 0.000003330622496

1e-012 0.000000000000000 0.000444049689291 0.000444049689291

1e-014 0.000000000000000 0.032949132854297 0.032949132854297

1e-016 0.000000000000000 1.000000000000000 1.000000000000000

>> A= [1 + 10^-14 2;1 2]

A =

1.0000 2.0000

>> format long

>> A= [1 + 10^-14 2;1 2]

A =

1.00000000000001 2.00000000000000

1.00000000000000 2.00000000000000

>> [Q,R] = cgs(A)

Q =

0.70710678118655 -0.68342428808113

0.70710678118654 0.73002139863212

R =

1.41421356237310 2.82842712474619

0 0.00000000000001

>> q1=Q(:,1)

q1 =

0.70710678118655

0.70710678118654

>> q2=Q(:,2)

q2 =

-0.68342428808113

0.73002139863212

>> q1'*q1

ans =

1.00000000000000

>> q2'*q2

ans =

1.00000000000000

>> q1'*q2

ans =

0.03294913285430

>> plot(q1(1),q1(2),'r*',q2(1),q2(2),'b*')

>> plot(q1(1),q1(2),'r*',q2(1),q2(2),'b*',0,0,'k*')

*nota : hemos añadido un punto en el (0,0)

Con el programa demo_cgs_mgs.m:

% Numerical comparison:

% Classical Gram-Schmidt and

% modified Gram-Schmidt

% First we generate a random matrix A

% with eigenvalues spaced by a factor 2

% between 2^-1 and 2^-80

[U,X]= qr(randn(80));

S = diag(2.^(-1:-1:-80));

A = U*S*U';

E = eig(A);

% We use now classical Gram-Schmid cgs(A)

% and modified Gram-Schmidt mgs(A) to compute Q and R

[Q1,R1]=cgs(A);

[Q2,R2]=mgs(A);

% We represent the diagonal R elements

% in logarithmic scale for both matrices R1 and R2.

x = 1:80;

y1 = log2(diag(R1));

y2 = log2(diag(R2));

y6 = log2(-sort(-abs(E)));

y3 = -23*ones(1,80);

y4 = -55*ones(1,80);

y5 = -x;

plot(x,y1,'bo',x,y2,'rx',x,y6,'g*',x,y3,':',x,y4,':',x,y5,'-');

xlabel('R(i,i) elements in logarithmic scale');

en matlab:

>> demo_cgs_mgs

>> 2^-55

ans =

2.775557561562891e-017

11-05-06

AX = B

A = Q * R = triángular superior

Q columnas (ó filas) son una base ortonormal, se llama ortogonal.

Q' * Q = Q * Q' = I

Q' Q RX = Q' B ; como Q' Q = I RX = Q' B

5. Singular value decomposition

(Página 7 del capitulo 5)

A = Q R = U S V'

Donde S es diagonal, y U y V' ortgonales

|s1 0 |

S = | s2 |

| 0 s3 |

Real Ortgonal

Compleja Unitaria

Ambas tmatrices tienen las columnas en bases ortonormales.

En marlab sería: >> [u,s,v] = svd(A) s1>>s2>>…>>sn

Hacemos en matlab:

>> help svd

SVD Singular value decomposition.

[U,S,V] = SVD(X) produces a diagonal matrix S, of the same

dimension as X and with nonnegative diagonal elements in

decreasing order, and unitary matrices U and V so that

X = U*S*V'.

S = SVD(X) returns a vector containing the singular values.

[U,S,V] = SVD(X,0) produces the "economy size"

decomposition. If X is m-by-n with m > n, then only the

first n columns of U are computed and S is n-by-n.

For m <= n, SVD(X,0) is equivalent to SVD(X).

[U,S,V] = SVD(X,'econ') also produces the "economy size"

decomposition. If X is m-by-n with m >= n, then it is

equivalent to SVD(X,0). For m < n, only the first m columns

of V are computed and S is m-by-m.

See also svds, gsvd.

Overloaded functions or methods (ones with the same name in other directories)

help sym/svd.m

Reference page in Help browser doc svd

>> A = rand(3)

A =

0.9501 0.4860 0.4565

0.2311 0.8913 0.0185

0.6068 0.7621 0.8214

>> [u,s,v]= svd(A)

u =

-0.6068 0.4443 -0.6591

-0.4007 -0.8871 -0.2290

-0.6865 0.1251 0.7163

s =

1.8109 0 0

0 0.6319 0

0 0 0.3748

v =

-0.5995 0.4637 -0.6523

-0.6489 -0.7587 0.0572

-0.4684 0.4576 0.7558

>> >> u*u'

ans =

1.0000 0 -0.0000

0 1.0000 0.0000

-0.0000 0.0000 1.0000

>> A - u*s*v'

ans =

1.0e-015 *

0.1110 -0.3886 0.1110

-0.3608 0 -0.3192

-0.3331 -0.6661 -0.1110

>> norm(A - u*s*v')

ans =

8.2930e-016

>> A = [1 0 0; 0 0 1; 1 1 1]

A =

1 0 0

0 0 1

1 1 1

>> spy(A)

>> A = [1 1 1 1; 1 1 1 1; 1 1 0 0; 1 1 1 1;1 1 0 0;1 1 1 1 ;1 1 1 1]

A =

1 1 1 1

1 1 0 0

1 1 1 1

1 1 0 0

1 1 1 1

>> spy(A)

En editor:

function L = letterL

L = [1 1 0 0 0 0;

1 1 0 0 0 0;

1 1 1 1 1 1;

1 1 1 1 1 1];

function O = letterO

O = [1 1 1 1 1 1;

1 1 1 1 1 1;

1 1 0 0 1 1;

1 1 1 1 1 1;

1 1 1 1 1 1];

function E = letterE

E = [1 1 1 1 1 1;

1 1 1 1 1 1;

1 1 1 0 0 0;

1 1 1 1 1 1;

1 1 1 0 0 0;

1 1 1 1 1 1 ;

1 1 1 1 1 1];

function H = letterH

H = [1 1 0 0 1 1;

1 1 0 0 1 1;

1 1 1 1 1 1;

1 1 0 0 1 1;

1 1 0 0 1 1];

Hacemos el programa:

%este script nos escribe la palabra "HELLO"

% por una matriz formada con 1 y 0

function h = hello

h = zeros(12,42);

h(3:10,3:8)= h(3:10,3:8) + letterH;

h(3:10,11:16)= h(3:10,11:16) + letterE;

h(3:10,19:24)= h(3:10,19:24) + letterL;

h(3:10,27:32)= h(3:10,27:32) + letterL;

h(3:10,35:40)= h(3:10,35:40) + letterO;

>> spy(hello)

>> [U,S,V]= svd(hello);

>> diag(S)

ans =

12.0087

3.7201

2.1501

1.8247

0.0000

>> pcolor(hello)

U1 =

0 0 0 0

-0.3598 -0.0869 0.1134 0.5917

-0.2923 -0.2906 -0.5635 -0.1124

-0.3483 -0.3132 0.4093 -0.3362

-0.2923 -0.2906 -0.5635 -0.1124

-0.4047 0.5567 -0.0461 -0.1555

0 0 0 0

>> V1 =V(:,1:4);

>> S1= S(1:4,1:4);

>> S1

S1 =

12.0087 0 0 0

0 3.7201 0 0

0 0 2.1501 0

0 0 0 1.8247

>> hello = U1*S1*V1';

>> pcolor(hello)

>> U2 = U(:,1:2);

>> V2 =V(:,1:2);

>> S2= S(1:2,1:2);

>> hello3 = U2*S2*V2'

16/05/06

(Tema 6)

CHAPTER 6. Quadratures

Tenemos el programa PesosNC:

% Este fichero almacena los

% pesos de las reglas de

% cuadratura de Newton-Cotes

% con m puntos

% m es un entero que satisface 2 <= m <= 11

% p es un vector columna que contiene

% los pesos para la regla de cuadratura

% de Newton-Cotes con m puntos.

function p = pesosNC(m);

if m==2

p=[1 1]'/2;

elseif m==3

p=[1 4 1]'/6;

elseif m==4

p=[1 3 3 1]'/8;

elseif m==5

p=[7 32 12 32 7]'/90;

elseif m==6

p=[19 75 50 50 75 19]'/288;

elseif m==7

p=[41 216 27 272 27 216 41]'/840;

elseif m==8

p=[751 3577 1323 2989 2989 1323 3577 751]'/17280;

elseif m==9

p=[989 5888 -928 10496 -4540 10496 -928 5888 989]'/28350;

elseif m==10

p=[2857 15741 1080 19344 5778 5778 19344 1080 15741 2857]'/89600;

else

p=[16067 106300 -48525 272400 -260550 427368 -260550 272400 -48525 106300 16067]'/598752;

end;

regla se simpson

en matlab:

>> pesosnc(3)

ans =

0.1667

0.6667

0.1667

>> pesosnc(5)

ans =

0.0778

0.3556

0.1333

0.3556

0.0778

Con el programa simpson:

% This function computes the

% Simpson quadrature rule for

% computing the integral of f(x) on

% the interval [a,b]

function I = simpson(fname,a,b)

fa = feval(fname,a);

fm = feval(fname, (a+b)/2);

fb = feval(fname, b);

I = ( fa + 4*fm + fb)*(b-a)/6;

Y el programa fname:

function y = fname(x)

y = 2 - 3*x + 4*x.^2;

En matlab:

>> I = simpson('fname',1,4)

I =

67.5000

Cambiamos el programa:

function y = fname(x)

%y = 2 - 3*x + 4*x.^2;

y = 1 - x.^4;

>> I = simpson('fname',0,2)

I =

-4.6667

>> 22/5

ans =

4.4000

Otro programa: simpsoncom:

% This function computes

% the Simpson composite quadrature rule

% with n subintervals

function I = simpsoncom(fname,a,b,n);

h = (b-a)/n;

x = [a:h:b];

f = feval(fname, x);

I = 0;

for k = 1:n

I = I + simpson('fname',x(k),x(k+1));

End

Y cuadraturaNC:

% Este fichero calcula la cuadratura

% de Newton-Cotes de la funcion fname

% con m puntos y extremos del intervalo a y b.

% a,b numeros reales

% m entero entre 2 <= m <=11

% El computo devuelve el numero I

function I = cuadraturaNC(fname,a,b,m)

p = pesosNC(m);

%x = linspace(a,b,m)';

h = (b-a)/(m-1);

x = (a:h:b)';

f = feval(fname,x);

I = (b-a)*(p'*f);

>> I = cuadraturaNC('fname',0,2,3)

Warning: Function call pesosNC invokes inexact match E:\computacion estadistica\16-05-06\pesosnc.m.

> In cuadraturaNC at 9

I =

-4.6667

>> I = cuadraturaNC('fname',0,2,4)

I =

-4.5185

>> I = cuadraturaNC('fname',0,2,5)

I =

-4.40000000000000

>> 22/5

ans =

4.40000000000000

Cambiamos fname:

function y = fname(x)

%y = 2 - 3*x + 4*x.^2;

%y = 1 - x.^4;

y = sin(x);

>> I = cuadraturaNC('fname',0,pi,3)

I =

2.09439510239320

>> I = cuadraturaNC('fname',0,pi,5)

I =

1.99857073182384

>> I = cuadraturaNC('fname',0,pi,11)

I =

2.00000000114677

>> I = simpsoncom('fname',0,pi,2)

I =

2.00455975498442

>> I = simpsoncom('fname',0,pi,10)

I =

2.00000678444180

>> I = simpsoncom('fname',0,pi,100)

I =

2.00000000067647

>> I = simpsoncom('fname',0,pi,1000)

I =

2.00000000000007

>> I = simpsoncom('fname',0,pi,10000)

I =

*nota, para borrar la pantalla: clear clc

Tenemos cuadraturacom;

% Este fichero calculo la cuadratura

% de Newton-Cotes compuesta con m nodos

% y n subintervalos

% a,b numeros reales, a < b

% m numero de puntos, 2 <= m <=11

% n numero de subintervalos

% El programa calcula

% I La aproximacion de Newton-Cotes compuesta

% de la integral de f(x) entre a y b.

% La regla es aplicada a cada uno de los

% n subintervalos iguales de [a,b].

function I = cuadcompNC(fname,a,b,m,n)

Delta = (b-a)/n;

h = Delta/(m-1);

x = a+h*(0:(n*(m-1)))';

p = pesosNC(m);

%x = linspace(a,b,n*(m-1)+1)';

f = feval(fname,x);

I = 0;

first = 1;

last = m;

for i=1:n

%A-adimos al producto interno el valor en

% el subintervalo i-esimo.

I = I + p'*f(first: last);

first = last;

last = last+m-1;

end

I = Delta*I;

Cambio fname:

function y = fname(x)

%y = 2 - 3*x + 4*x.^2;

%y = 1 - x.^4;

%y = sin(x);

y = sin(x)./x;

con el programa cuadraturacomNC:

% Este fichero calcula la cuadratura

% de Newton-Cotes de la funcion fname

% con m puntos y extremos del intervalo a y b.

% a,b numeros reales

% m entero entre 2 <= m <=11

% El computo devuelve el numero I

function I = cuadraturaNC(fname,a,b,m)

p = pesosNC(m);

%x = linspace(a,b,m)';

h = (b-a)/(m-1);

x = (a:h:b)';

f = feval(fname,x);

I = (b-a)*(p'*f);

18/05/06

Tema 6 del book moler

Ejercicio 5.9

5.9. Statistical Reference Datasets. NIST, the National Institute of Standards and

Technology, is the branch of the U.S. Department of Commerce responsible

for setting national and international standards. NIST maintains Statistical

Reference Datasets, StRD, for use in testing and certifying statistical soft-

ware. The home page on the Web is [3]. Data sets for linear least squares are

under \Linear Regression." This exercise involves two of the NIST reference

data sets.

² Norris. Linear polynomial for calibration of ozone monitors.

² Pontius. Quadratic polynomial for calibration of load cells.

For each of these data sets, follow the Web links labeled

² Data ¯le (ASCII Format)

² Certi¯ed Values

² Graphics

Download each ASCII ¯le. Extract the observations. Compute the polyno-

mial coe±cients. Compare the coe±cients with the certi¯ed values. Make

plots similar to the NIST plots of both the ¯t and the residuals.

Copiamos de: http://www.itl.nist.gov/div898/strd/lls/data/LINKS/DATA/Norris.dat

NIST/ITL StRD

Dataset Name: Norris (Norris.dat)

File Format: ASCII

Certified Values (lines 31 to 46)

Data (lines 61 to 96)

Procedure: Linear Least Squares Regression

Reference: Norris, J., NIST.

Calibration of Ozone Monitors.

Data: 1 Response Variable (y)

1 Predictor Variable (x)

36 Observations

Lower Level of Difficulty

Observed Data

Model: Linear Class

2 Parameters (B0,B1)

y = B0 + B1*x + e

Certified Regression Statistics

Standard Deviation

Parameter Estimate of Estimate

B0 -0.262323073774029 0.232818234301152

B1 1.00211681802045 0.429796848199937E-03

Residual

Standard Deviation 0.884796396144373

R-Squared 0.999993745883712

Certified Analysis of Variance Table

Source of Degrees of Sums of Mean

Variation Freedom Squares Squares F Statistic

Regression 1 4255954.13232369 4255954.13232369 5436385.54079785

Residual 34 26.6173985294224 0.782864662630069

Data: y x

0.1 0.2

338.8 337.4

118.1 118.2

888.0 884.6

9.2 10.1

228.1 226.5

668.5 666.3

998.5 996.3

449.1 448.6

778.9 777.0

559.2 558.2

0.3 0.4

0.1 0.6

778.1 775.5

668.8 666.9

339.3 338.0

448.9 447.5

10.8 11.6

557.7 556.0

228.3 228.1

998.0 995.8

888.8 887.6

119.6 120.2

0.3 0.3

0.6 0.3

557.6 556.8

339.3 339.1

888.0 887.2

998.5 999.0

778.9 779.0

10.2 11.1

117.6 118.3

228.9 229.2

668.4 669.1

449.2 448.9

0.2 0.5

*Para importer datos: File Import Data seleccionas los datos que quieres importar

aceptas a todo y finalizas luego en matlab pones el nombre del archivo y se abre en el programa.

En matlab:

>> y = norris2(:,1)

y =

0.1000

338.8000

118.1000

888.0000

9.2000

228.1000

668.5000

998.5000

449.1000

778.9000

559.2000

0.3000

0.1000

778.1000

668.8000

339.3000

448.9000

10.8000

557.7000

228.3000

998.0000

888.8000

119.6000

0.3000

0.6000

557.6000

339.3000

888.0000

998.5000

778.9000

10.2000

117.6000

228.9000

668.4000

449.2000

0.2000

>> x = norris2(:,2)

x =

0.2000

337.4000

118.2000

884.6000

10.1000

226.5000

666.3000

996.3000

448.6000

777.0000

558.2000

0.4000

0.6000

775.5000

666.9000

338.0000

447.5000

11.6000

556.0000

228.1000

995.8000

887.6000

120.2000

0.3000

556.8000

339.1000

887.2000

999.0000

779.0000

11.1000

118.3000

229.2000

669.1000

448.9000

0.5000

>> n = length(x)

n =

>> A = ones(36,2);

>> A(:,2) = x;

>> A

A =

1.0000 0.2000

1.0000 337.4000

1.0000 118.2000

1.0000 884.6000

1.0000 10.1000

1.0000 226.5000

1.0000 666.3000

1.0000 996.3000

1.0000 448.6000

1.0000 777.0000

1.0000 558.2000

1.0000 0.4000

1.0000 0.6000

1.0000 775.5000

1.0000 666.9000

1.0000 338.0000

1.0000 447.5000

1.0000 11.6000

1.0000 556.0000

1.0000 228.1000

1.0000 995.8000

1.0000 887.6000

1.0000 120.2000

1.0000 0.3000

1.0000 556.8000

1.0000 339.1000

1.0000 887.2000

1.0000 999.0000

1.0000 779.0000

1.0000 11.1000

1.0000 118.3000

1.0000 229.2000

1.0000 669.1000

1.0000 448.9000

1.0000 0.5000

>> rank(A)

ans =

>> rank([A y])

ans =

>> beta = A\y

beta =

-0.2623

1.0021

>> format long

>> beta = A\y

beta =

-0.26232307377407

1.00211681802045

Ahora pegamos de la misma página :

Pointius:

NIST/ITL StRD

Dataset Name: Pontius

File Format: ASCII

Certified Values (lines 31 to 47)

Data (lines 61 to 100)

Procedure: Linear Least Squares Regression

Reference: Pontius, P., NIST.

Load Cell Calibration.

Data: 1 Response Variable (y)

1 Predictor Variable (x)

40 Observations

Lower Level of Difficulty

Observed Data

Model: Quadratic Class

3 Parameters (B0,B1,B2)

y = B0 + B1*x + B2*(x**2)

Certified Regression Statistics

Standard Deviation

Parameter Estimate of Estimate

B0 0.673565789473684E-03 0.107938612033077E-03

B1 0.732059160401003E-06 0.157817399981659E-09

B2 -0.316081871345029E-14 0.486652849992036E-16

Residual

Standard Deviation 0.205177424076185E-03

R-Squared 0.999999900178537

Certified Analysis of Variance Table

Source of Degrees of Sums of Mean

Variation Freedom Squares Squares F Statistic

Regression 2 15.6040343244198 7.80201716220991 185330865.995752

Residual 37 0.155761768796992E-05 0.420977753505385E-07

Data: y x

.11019 150000

.21956 300000

.32949 450000

.43899 600000

.54803 750000

.65694 900000

.76562 1050000

.87487 1200000

.98292 1350000

1.09146 1500000

1.20001 1650000

1.30822 1800000

1.41599 1950000

1.52399 2100000

1.63194 2250000

1.73947 2400000

1.84646 2550000

1.95392 2700000

2.06128 2850000

2.16844 3000000

.11052 150000

.22018 300000

.32939 450000

.43886 600000

.54798 750000

.65739 900000

.76596 1050000

.87474 1200000

.98300 1350000

1.09150 1500000

1.20004 1650000

1.30818 1800000

1.41613 1950000

1.52408 2100000

1.63159 2250000

1.73965 2400000

1.84696 2550000

1.95445 2700000

2.06177 2850000

2.16829 3000000

Y en matlab se importan igual los datos y (solo los numeros le llamo pointius2):

pointius2 =

1.0e+006 *

0.00000011019000 0.15000000000000

0.00000021956000 0.30000000000000

0.00000032949000 0.45000000000000

0.00000043899000 0.60000000000000

0.00000054803000 0.75000000000000

0.00000065694000 0.90000000000000

0.00000076562000 1.05000000000000

0.00000087487000 1.20000000000000

0.00000098292000 1.35000000000000

0.00000109146000 1.50000000000000

0.00000120001000 1.65000000000000

0.00000130822000 1.80000000000000

0.00000141599000 1.95000000000000

0.00000152399000 2.10000000000000

0.00000163194000 2.25000000000000

0.00000173947000 2.40000000000000

0.00000184646000 2.55000000000000

0.00000195392000 2.70000000000000

0.00000206128000 2.85000000000000

0.00000216844000 3.00000000000000

0.00000011052000 0.15000000000000

0.00000022018000 0.30000000000000

0.00000032939000 0.45000000000000

0.00000043886000 0.60000000000000

0.00000054798000 0.75000000000000

0.00000065739000 0.90000000000000

0.00000076596000 1.05000000000000

0.00000087474000 1.20000000000000

0.00000098300000 1.35000000000000

0.00000109150000 1.50000000000000

0.00000120004000 1.65000000000000

0.00000130818000 1.80000000000000

0.00000141613000 1.95000000000000

0.00000152408000 2.10000000000000

0.00000163159000 2.25000000000000

0.00000173965000 2.40000000000000

0.00000184696000 2.55000000000000

0.00000195445000 2.70000000000000

0.00000206177000 2.85000000000000

0.00000216829000 3.00000000000000

Nota : si en el programa no tenemos la opción de “import data” , importaríamkos los datos de la siguiente manera: Pointius = load(`pointius.txt');

>> y = pointius2(:,1);

>> x = pointius2(:,2);

>> n = lenght(x)

>> n = length(x)

n =

>> A = ones(40,3);

>> A(:,2) = x;

>> A(:,3) = x.^2;

A =

1.0e+012 *

0.00000000000100 0.00000015000000 0.02250000000000

0.00000000000100 0.00000030000000 0.09000000000000

0.00000000000100 0.00000045000000 0.20250000000000

0.00000000000100 0.00000060000000 0.36000000000000

0.00000000000100 0.00000075000000 0.56250000000000

0.00000000000100 0.00000090000000 0.81000000000000

0.00000000000100 0.00000105000000 1.10250000000000

0.00000000000100 0.00000120000000 1.44000000000000

0.00000000000100 0.00000135000000 1.82250000000000

0.00000000000100 0.00000150000000 2.25000000000000

0.00000000000100 0.00000165000000 2.72250000000000

0.00000000000100 0.00000180000000 3.24000000000000

0.00000000000100 0.00000195000000 3.80250000000000

0.00000000000100 0.00000210000000 4.41000000000000

0.00000000000100 0.00000225000000 5.06250000000000

0.00000000000100 0.00000240000000 5.76000000000000

0.00000000000100 0.00000255000000 6.50250000000000

0.00000000000100 0.00000270000000 7.29000000000000

0.00000000000100 0.00000285000000 8.12250000000000

0.00000000000100 0.00000300000000 9.00000000000000

0.00000000000100 0.00000015000000 0.02250000000000

0.00000000000100 0.00000030000000 0.09000000000000

0.00000000000100 0.00000045000000 0.20250000000000

0.00000000000100 0.00000060000000 0.36000000000000

0.00000000000100 0.00000075000000 0.56250000000000

0.00000000000100 0.00000090000000 0.81000000000000

0.00000000000100 0.00000105000000 1.10250000000000

0.00000000000100 0.00000120000000 1.44000000000000

0.00000000000100 0.00000135000000 1.82250000000000

0.00000000000100 0.00000150000000 2.25000000000000

0.00000000000100 0.00000165000000 2.72250000000000

0.00000000000100 0.00000180000000 3.24000000000000

0.00000000000100 0.00000195000000 3.80250000000000

0.00000000000100 0.00000210000000 4.41000000000000

0.00000000000100 0.00000225000000 5.06250000000000

0.00000000000100 0.00000240000000 5.76000000000000

0.00000000000100 0.00000255000000 6.50250000000000

0.00000000000100 0.00000270000000 7.29000000000000

0.00000000000100 0.00000285000000 8.12250000000000

0.00000000000100 0.00000300000000 9.00000000000000

>> rank(A)

ans =

>> rank([A y])

ans =

>> beta = A\y

beta =

1.0e-003 *

0.67356578947336

0.00073205916040

-0.00000000000316

1) Ecuaciones normales AX = B ; A' A X = A' B

A' A X = M; MX = C

2) A = Q R

RX = Q' R X = R(Q*__); donde (Q*___) = D

3) A = U S V'

X = VS + U' B

>> M = A'*A

M =

1.0e+026 *

0.00000000000000 0.00000000000000 0.00000000000129

0.00000000000000 0.00000000000129 0.00000297675000

0.00000000000129 0.00000297675000 7.31699325000000

>> C = A' * y

C =

1.0e+014 *

0.00000000000046

0.00000093646979

2.15689932675000

>> beta1= M\C

Warning: Matrix is close to singular or badly scaled.

Results may be inaccurate. RCOND = 4.938239e-027.

beta1 =

1.0e-003 *

0.67356578947410

0.00073205916040

-0.00000000000316

>> [Q,R]=qr(A);

>> D = Q'*y;

>> beta2= R\D;

>> beta2

beta2 =

1.0e-003 *

0.67356578947465

0.00073205916040

-0.00000000000316

>> [U,S,V]= svd(A);

>> V

V =

0.00000000000018 0.00000129952456 0.99999999999916

0.00000040682694 0.99999999999907 -0.00000129952456

0.99999999999992 -0.00000040682694 0.00000000000035

>> S1 = pinv(S);

>> S1(1,2)

ans =

>> diag(S1)

ans =

0.00000000000004

0.00000035250209

0.52607450611596

>> beta3 = V*S1*U'*y;

>> beta3

beta3 =

1.0e-003 *

0.67356576935354

0.00073205916049

-0.00000000000316

>> norm(A*(beta2 - beta3))

ans =

2.128137835323598e-009

>> norm(A'*A*(beta2 - beta3))

ans =

5.735906807875822e+004

Cogemos el ejercicio 5.10 del moler ( y es el que entregamos de práctica)

5.10. Filip data set. One of the Statistical Reference Datasets from the NIST is the

\Filip" dataset. The data consists of several dozen observations of a variable

y at di®erent values of x. The task is to model y by a polynomial of degree

10 in x.

This dataset is controversial. A search of the Web for \¯lip strd" will ¯nd

several dozen postings, including the original page at NIST [3]. Some math-

ematical and statistical packages are able to reproduce the polynomial coef-

¯cients that NIST has decreed to be the \certi¯ed values." Other packages

give warning or error messages that the problem is too badly conditioned to

solve. A few packages give di®erent coe±cients without warning. The Web

o®ers several opinions about whether or not this is a reasonable problem.

Let's see what MATLAB does with it.

The data set is available from the NIST Web site. There is one line for

each data point. The data is given with the ¯rst number on the line a value

of y, and the second number the corresponding x. The x-values are not

monotonically ordered, but it is not necessary to sort them. Let n be the

number of data points and p = 11 the number of polynomial coe±cients.

(a) As your ¯rst experiment, load the data into MATLAB, plot it with '.'

as the line type, and then invoke the Basic Fitting tool available under the

Tools menu on the ¯gure window. Select the 10th degree polynomial ¯t. You

will be warned that the polynomial is badly conditioned, but ignore that for

now. How do the computed coe±cients compare with the certi¯ed values on

the NIST Web page? How does the plotted ¯t compare with the graphic on

the NIST Web page? The basic ¯tting tool also displays the norm of the

residuals, krk. Compare this with the NIST quantity \Residual Standard

Deviation," which is

krk pn ¡ p

(b) Examine this data set more carefully by using six di®erent methods to

compute the polynomial ¯t. Explain all the warning messages you receive

during these computations.

² Poly¯t: Use polyfit(x,y,10)

² Backslash: Use Xny where X is the n-by-p truncated Vandermonde ma-

trix with elements

Xi;j = xp¡j

i ; i = 1; : : : ; n; j = 1; : : : ; p

² Pseudoinverse: Use pinv(X)*y

² Normal equations: Use inv(X'*X)*X'*y.

² Centering: Let ¹ = mean(x); ¾ = std(x); t = (x ¡ ¹)=¾.

Use polyfit(t,y,10).

² Certi¯ed coe±cients: Obtain the coe±cients from the NIST Web page.

di®erent methods?

24 Chapter 5. Least Squares

(d) Which one of the six methods gives a very poor ¯t? (Perhaps the packages

that are criticized on theWeb for reporting bad results are using this method.)

(e) Plot the ¯ve good ¯ts. Use dots, '.', at the data values and curves

obtained by evaluating the polynomials at a few hundred points over the

range of the x's. The plot should look like ¯gure 5.5. There are ¯ve di®erent

plots, but only two visually distinct ones. Which methods produce which

plots?

(f) Why do poly¯t and backslash give di®erent results?

"9 "8 "7 "6 "5 "4 "3

0.76

0.78

0.8

0.82

0.84

0.86

0.88

0.9

0.92

0.94

NIST Filip data set

Data

Rank 11

Rank 10

Figure 5.5. NIST Filip standard reference data set

TEMA 6 del mooler:

Ejercicio 6.4:

6.4. Use quadtx with various tolerances to compute pi by approximating

pi = integral desde -1 a 1 de : 2/(1+x^2) dx

solución : 2*arcotg entre -1 y 1 = 2 [ pi/4 - (-pi/4)] = pi

How do the accuracy and the function evaluation count vary with tolerance?

Una manera de calcular pi es por cuadraturas, pero no es muy eficiente.

F(x) = 2/(1+x^2)

Fname.m

Inline

En matlab:

>> F = inline('2./(1+x^2)')

F =

Inline function:

F(x) = 2./(1+x^2)

>> G = inline('x.^y','x','y')

G =

Inline function:

G(x,y) = x.^y

Con el programa quadtx:

function [Q,fcount] = quadtx(F,a,b,tol,varargin)

%QUADTX Evaluate definite integral numerically.

% Q = QUADTX(F,A,B) approximates the integral of F(x) from A to B

% to within a tolerance of 1.e-6. F is a string defining a function

% of a single variable, an inline function, a function handle, or a

% symbolic expression involving a single variable.

% Q = QUADTX(F,A,B,tol) uses the given tolerance instead of 1.e-6.

% Arguments beyond the first four, Q = QUADTX(F,a,b,tol,p1,p2,...),

% are passed on to the integrand, F(x,p1,p2,..).

% [Q,fcount] = QUADTX(F,...) also counts the number of evaluations

% of F(x).

% See also QUAD, QUADL, DBLQUAD, QUADGUI.

% Make F callable by feval.

if ischar(F) & exist(F)~=2

F = inline(F);

elseif isa(F,'sym')

F = inline(char(F));

end

% Default tolerance

if nargin < 4 | isempty(tol)

tol = 1.e-6;

end

% Initialization

c = (a + b)/2;

fa = feval(F,a,varargin{:});

fc = feval(F,c,varargin{:});

fb = feval(F,b,varargin{:});

% Recursive call

[Q,k] = quadtxstep(F, a, b, tol, fa, fc, fb, varargin{:});

fcount = k + 3;

% ---------------------------------------------------------

function [Q,fcount] = quadtxstep(F,a,b,tol,fa,fc,fb,varargin)

% Recursive subfunction used by quadtx.

h = b - a;

c = (a + b)/2;

fd = feval(F,(a+c)/2,varargin{:});

fe = feval(F,(c+b)/2,varargin{:});

Q1 = h/6 * (fa + 4*fc + fb);

Q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);

if abs(Q2 - Q1) <= tol

Q = Q2 + (Q2 - Q1)/15;

fcount = 2;

else

[Qa,ka] = quadtxstep(F, a, c, tol, fa, fd, fc, varargin{:});

[Qb,kb] = quadtxstep(F, c, b, tol, fc, fe, fb, varargin{:});

Q = Qa + Qb;

fcount = ka + kb + 2;

end

>> I = quadtx(F,-1,1,10^-6)

I =

3.14159265370804

I2 = cuadcompNC(F,-1,1,3,10)

I2=

3.14159261393922

I2= cuadcmpNC(F,-1,1,8,10)

I2 =

3.14159265358986

I2= cuadcompNC(F,-1,1,8,100)

I2 =

3.14159265358979

Otro ejercicio es el 6.6 (página14 del capitulo 6 del mooler)

6.6. The error function, erf(x), is de¯ned by an integral,

erf(x) = 2/raiz(pi) * integral de 0 a x de “e” elevado a “-x^2”

Use quadtx to tabulate erf(x) for x = 0:1; 0:2; : : : ; 1:0: Compare the results

with the built-in Matlab function erf(x).

*Nota para el 6.8 . Tenemos una gamma(n+1) = n!

30-05-06

Capitulo 8 (chapter 7)

*Metodo de Picard

Tenemos el archivo:

Harmonic:

function ydot = harmonic(t,y)

ydot = [y(2); -y(1)];

en matlab:

>> [t,y] = ode23('harmonic',[0 5],[3;0]);

>> [t,y] = ode23('harmonic',[0 5],[3;0])

t =

0.0000

0.0002

0.0008

0.0042

0.0208

0.1008

0.2361

0.4164

0.6367

0.8973

1.1880

1.4396

1.5862

1.7328

1.8913

2.0918

2.3317

2.6147

2.8857

3.0680

3.2062

3.3445

3.5157

3.7276

3.9793

4.2768

4.5361

4.6975

4.8589

5.0000

y =

3.0000 0

3.0000 -0.0001

3.0000 -0.0005

3.0000 -0.0025

3.0000 -0.0125

2.9993 -0.0625

2.9848 -0.3019

2.9167 -0.7019

2.7435 -1.2134

2.4118 -1.7834

1.8704 -2.3443

1.1194 -2.7813

0.3917 -2.9719

-0.0467 -2.9972

-0.4841 -2.9581

-0.9448 -2.8446

-1.4923 -2.5993

-2.0671 -2.1698

-2.5902 -1.5058

-2.8981 -0.7573

-2.9872 -0.2194

-2.9889 0.1945

-2.9335 0.6046

-2.7875 1.0955

-2.4945 1.6571

-2.0029 2.2258

-1.2619 2.7144

-0.5234 2.9467

-0.0431 2.9924

0.4383 2.9604

0.8503 2.8692

>> plot(t,y(:,1))

> [t,y] = ode23('harmonic',[0 50],[3;0]);

>> plot(t,y(:,1))

Tenemos el programa :

Harmonic_damped al que llamamos por comodidad : harmonicd :

function ydot = harmonic_damped(t,y)

alpha = .5;

ydot = [y(2); -y(1) - alpha *y(2)];

en matlab:

>> [t,y] = ode23('harmonicd',[0 50],[3;0]);

>> plot(t,y(:,1))

Cambiamos el alfa del programa:

function ydot = harmonic_damped(t,y)

alpha = .1;

ydot = [y(2); -y(1) - alpha *y(2)];

>> [t,y] = ode23('harmonicd',[0 50],[3;0]);

>> plot(t,y(:,1))

Cambiamos el programa:

function ydot = harmonic_damped(t,y)

alpha = .1; beta = 1.2;

ydot = [y(2); -y(1) - alpha *(y(2).^beta)];

en matlab:

>> [t,y] = ode23('harmonicd',[0 50],[3;0]);

>> plot(t,y(:,1))

Warning: Imaginary parts of complex X and/or Y arguments ignored.

>> plot(y(:,1),y(:,2))

Warning: Imaginary parts of complex X and/or Y arguments ignored.

>> [t,y] = ode23('harmonicd',[0 50],[7;2]);

>> [t,y] = ode23('harmonicd',[0 100],[7;2]);

>> plot(y(:,1),y(:,2))

Warning: Imaginary parts of complex X and/or Y arguments ignored.

Tenemos el programa:

Twobody:

function ydot = twobody(t,y)

r = sqrt(y(1)^2 + y(2)^2);

ydot = [y(3); y(4); -y(1)/r^3 ; -y(2)/r^3];

>> [t,y] = ode23('twobody',[0 3],[5;0;0;0.1]);

>> plot(y(:,1),y(:,2))

>> [t,y] = ode23('twobody',[0 10],[10;0;0;0.1]);

>> plot(y(:,1),y(:,2))

>> [t,y] = ode23('twobody',[0 10],[3;0;0;0.1]);

>> plot(y(:,1),y(:,2))

>> [t,y] = ode23('twobody',[0 20],[3;0;0;0.1]);

>> plot(y(:,1),y(:,2))

>> [t,y] = ode45('twobody',[0 20],[3;0;0;0.1]);

>> plot(y(:,1),y(:,2))

>> [t,y] = ode45('twobody',[0 200],[3;0;0;0.1]);

>> plot(y(:,1),y(:,2))

1-junio-06

Ejercicio del libro moler (26 Chapter 5. Least Squares):

5.12. Planetary orbit [2]. The expression z = ax2+bxy+cy2+dx+ey+f is known

as a quadratic form. The set of points (x; y) where z = 0 is a conic section.

It can be an ellipse, a parabola, or a hyperbola, depending on the sign of

the discriminant b2 ¡ 4ac. Circles and lines are special cases. The equation

z = 0 can be normalized by dividing the quadratic form by any nonzero

coe±cient. For example, if f 6= 0, we can divide all the other coe±cients by

f and obtain a quadratic form with the constant term equal to one. You can

use the Matlab meshgrid and contour functions to plot conic sections. Use

meshgrid to create arrays X and Y. Evaluate the quadratic form to produce

Z. Then use contour to plot the set of points where Z is zero.

[X,Y] = meshgrid(xmin:deltax:xmax,ymin:deltay:ymax);

Z = a*X.^2 + b*X.*Y + c*Y.^2 + d*X + e*Y + f;

contour(X,Y,Z,[0 0])

A planet follows an elliptical orbit. Here are ten observations of its position

in the (x; y) plane:

x = [1.02 .95 .87 .77 .67 .56 .44 .30 .16 .01]';

y = [0.39 .32 .27 .22 .18 .15 .13 .12 .13 .15]';

(a) Determine the coe±cients in the quadratic form that ¯ts this data in

the least squares sense by setting one of the coe±cients equal to one and

solving a 10-by-5 overdetermined system of linear equations for the other

¯ve coe±cients. Plot the orbit with x on the x-axis and y on the y-axis.

Superimpose the ten data points on the plot.

(b) This least squares problem is nearly rank de¯cient. To see what e®ect this

has on the solution, perturb the data slightly by adding to each coordinate

of each data point a random number uniformly distributed in the interval

[-.005,.005]. Compute the new coe±cients resulting from the perturbed data.

Plot the new orbit on the same plot with the old orbit. Comment on your

comparison of the sets of coe±cients and the orbits.

En matlab:

Practica3

[x,y]=meshgrid(-2:0.01:2,-2:0.01:2);

a=1;b=1;c=1;e=1;f=-1;

>> z = a*x.^2+b*x.*y+c*y.^2+b*x+e*y+f;

>> contour(x,y,z,[0 0])

>> contour(x,y,z,5)

>> contour(x,y,z,10)

'Programación con Matlab'

Para cada punto de x e y tenemos una ecuación diferente de:

= ax^2 + bxy +cy^2+dx+ey

C= (a b c d e)' ; A*C = (1 …. 1)' = r

X`2= (x1^2 …. X10^2)'

A=[X`2 XY Y^2 X Y]

C = A/r

En matlab:

>> x = [1.02 .95 .87 .77 .67 .56 .44 .30 .16 .01]';

>> y = [0.39 .32 .27 .22 .18 .15 .13 .12 .13 .15]';

>> x

x =

1.0200

0.9500

0.8700

0.7700

0.6700

0.5600

0.4400

0.3000

0.1600

0.0100

>> y

y =

0.3900

0.3200

0.2700

0.2200

0.1800

0.1500

0.1300

0.1200

0.1300

0.1500

>> A = [x.^2 x.*y y.^2 x y];

>> size(A)

ans =

10 5

>> r= ones(10,1)

r =

>> [X,Y]= meshgrid(-2:0.01:2,-2:0.01:2);

>> z = a*X.^2+b*X.*Y+c*Y.^2+b*X+e*Y+f;

>> A = [x.^2 x.*y y.^2 x y];

>> C=A\r;

>> a = C(1);b=C(2);c= C(3);d=C(4);e=C(5);

>> Z = a*X.^2+b*X.*Y+c*Y.^2+b*X+e*Y+f;

>> [X,Y]= meshgrid(-2:0.01:2,-2:0.01:2);

>> x = [1.02 .95 .87 .77 .67 .56 .44 .30 .16 .01]';

>> y = [0.39 .32 .27 .22 .18 .15 .13 .12 .13 .15]';

>> Z = a*X.^2+b*X.*Y+c*Y.^2+b*X+e*Y+f;

>> C=A\r;

>> a = C(1);b=C(2);c= C(3);d=C(4);e=C(5);

>> contour(X,Y,Z,[0 0])

>> plot(x,y,'r*')

>> x = x + 0.005*(rand(10,1)-0.5);

>> y = y + 0.005*(rand(10,1)-0.5);

>> A = [x.^2 x.*y y.^2 x y];

>> C = A\r;

>> Z = a*X.^2+b*X.*Y+c*Y.^2+b*X+e*Y+f;

>> a = C(1); b = C(2); c = C(3); d = C(4); e = C(5);

>> Z = a*X.^2+b*X.*Y+c*Y.^2+b*X+e*Y+f;

>> contour(X,Y,Z,[0 0])

>> hold on

PROGRAMAS USADOS:

BSLASHTX

function x = bslashtx(A,b)

% BSLASHTX Solve linear system (backslash)

% x = bslashtx(A,b) solves A*x = b

[n,n] = size(A);

if isequal(triu(A,1),zeros(n,n))

% Lower triangular

x = forward(A,b);

return

elseif isequal(tril(A,-1),zeros(n,n))

% Upper triangular

x = backsubs(A,b);

return

elseif isequal(A,A')

[R,fail] = chol(A);

if ~fail

% Positive definite

y = forward(R',b);

x = backsubs(R,y);

return

end

% Triangular factorization

[L,U,p] = lutx(A);

% Permutation and forward elimination

y = forward(L,b(p));

% Back substitution

x = backsubs(U,y);

% ------------------------------

function x = forward(L,x)

% FORWARD. Forward elimination.

% For lower triangular L, x = forward(L,b) solves L*x = b.

[n,n] = size(L);

x(1) = x(1)/L(1,1);

for k = 2:n

j = 1:k-1;

x(k) = (x(k) - L(k,j)*x(j))/L(k,k);

end

% ------------------------------

function x = backsubs(U,x)

% BACKSUBS. Back substitution.

% For upper triangular U, x = backsubs(U,b) solves U*x = b.

[n,n] = size(U);

x(n) = x(n)/U(n,n);

for k = n-1:-1:1

j = k+1:n;

x(k) = (x(k) - U(k,j)*x(j))/U(k,k);

end

CGS

% Gram-Schmidt orthogonalization (unstable)

function [Q,R] = cgs(A)

n = length(A);

% starting variables

R = zeros(n,n);

Q = zeros(n,n);

V = zeros(n,n);

% Normalizing the first column vector A

R(1,1)=norm(A(:,1));

Q(:,1)= A(:,1)/R(1,1);

% Classical Gram-Schmidt

for j = 2:n

V(:,j) = A(:,j);

for i = 1:j-1

R(i,j) = Q(:,i)'*A(:,j);

V(:,j)= V(:,j)-R(i,j)*Q(:,i);

end

R(j,j)=norm(V(:,j));

Q(:,j)=V(:,j)/R(j,j);

end

CUADCOMNC

% Este fichero calcula la cuadratura

% de Newton-Cotes de la funcion fname

% con m puntos y extremos del intervalo a y b.

% a,b numeros reales

% m entero entre 2 <= m <=11

% El computo devuelve el numero I

function I = cuadraturaNC(fname,a,b,m)

p = pesosNC(m);

%x = linspace(a,b,m)';

h = (b-a)/(m-1);

x = (a:h:b)';

f = feval(fname,x);

I = (b-a)*(p'*f);

CUADRATURA

% Este fichero calcula la cuadratura

% de Newton-Cotes de la funcion fname

% con m puntos y extremos del intervalo a y b.

% a,b numeros reales

% m entero entre 2 <= m <=11

% El computo devuelve el numero I

function I = cuadraturaNC(fname,a,b,m)

p = pesosNC(m);

%x = linspace(a,b,m)';

h = (b-a)/(m-1);

x = (a:h:b)';

f = feval(fname,x);

I = (b-a)*(p'*f);

CUADRATURANC

% Este fichero calculo la cuadratura

% de Newton-Cotes compuesta con m nodos

% y n subintervalos

% a,b numeros reales, a < b

% m numero de puntos, 2 <= m <=11

% n numero de subintervalos

% El programa calcula

% I La aproximacion de Newton-Cotes compuesta

% de la integral de f(x) entre a y b.

% La regla es aplicada a cada uno de los

% n subintervalos iguales de [a,b].

function I = cuadcompNC(fname,a,b,m,n)

Delta = (b-a)/n;

h = Delta/(m-1);

x = a+h*(0:(n*(m-1)))';

p = pesosNC(m);

%x = linspace(a,b,n*(m-1)+1)';

f = feval(fname,x);

I = 0;

first = 1;

last = m;

for i=1:n

%A-adimos al producto interno el valor en

% el subintervalo i-esimo.

I = I + p'*f(first: last);

first = last;

last = last+m-1;

end

I = Delta*I;

DEMO_CGS_MGS

% Numerical comparison:

% Classical Gram-Schmidt and

% modified Gram-Schmidt

% First we generate a random matrix A

% with eigenvalues spaced by a factor 2

% between 2^-1 and 2^-80

[U,X]= qr(randn(80));

S = diag(2.^(-1:-1:-80));

A = U*S*U';

E = eig(A);

% We use now classical Gram-Schmid cgs(A)

% and modified Gram-Schmidt mgs(A) to compute Q and R

[Q1,R1]=cgs(A);

[Q2,R2]=mgs(A);

% We represent the diagonal R elements

% in logarithmic scale for both matrices R1 and R2.

x = 1:80;

y1 = log2(diag(R1));

y2 = log2(diag(R2));

y6 = log2(-sort(-abs(E)));

y3 = -23*ones(1,80);

y4 = -55*ones(1,80);

y5 = -x;

plot(x,y1,'bo',x,y2,'rx',x,y6,'g*',x,y3,':',x,y4,':',x,y5,'-');

xlabel('R(i,i) elements in logarithmic scale');

DEMO_INTER

% This script shows the graphic of cos(4x)

% by interpolating with a polynomial of degree

% 12 found by least squares.

close all

t = (linspace(0,5,50))';

A = vander(t);

A = A(:,38:50);

b = (cos(4*t));

u = 0:0.01:5;

v = cos(4*u);

p = A\b;

pvals = horner(p,t);

plot(t,b,'or',t,pvals,'b',u,v,'g')

DEMO_ORTOG

% This script explores the lost of orthogonality

% due to the intrinsic numerical unstability

% of the classical and

% modified Gram-Schmidt algorithms.

clc

disp(' Losing orthogonality on Gram-Schmidt algorithms ')

disp(' delta Error on QR Error on CGS Error on MGS ')

disp('------------------------------------------------------------------------')

for delta = logspace(-4,-16,7)

A = [1+delta 2;1 2];

[Q1,R1]= qr(A);

[Q2,R2]= cgs(A);

[Q3,R3]= mgs(A);

error1 = norm(Q1'*Q1 - eye(2));

error2 = norm(Q2'*Q2 - eye(2));

error3 = norm(Q3'*Q3 - eye(2));

disp(sprintf(' %5.0e %20.15f %20.15f %20.15f ', delta,error1,error2,error3))

end

EULER

function y = euler(f,t0,y0,h,tf)

t = t0;

y = y0;

while t <= tf

y = y + h*feval(f,t,y);

t = t + h;

end

FZEROTX

function b = fzerotx(F,ab,varargin)

%FZEROTX Textbook version of FZERO.

% x = fzerotx(F,[a,b]) tries to find a zero of F(x) between a and b.

% F(a) and F(b) must have opposite signs. fzerotx returns one

% end point of a small subinterval of [a,b] where F changes sign.

% Arguments beyond the first two, fzerotx(F,[a,b],p1,p2,...),

% are passed on, F(x,p1,p2,..).

% Example:

% fzerotx('sin(x)',[1,4])

% Make F callable by feval.

if ischar(F) & exist(F)~=2

F = inline(F);

elseif isa(F,'sym')

F = inline(char(F));

end

% Initialize.

a = ab(1);

b = ab(2);

fa = feval(F,a,varargin{:});

fb = feval(F,b,varargin{:});

if sign(fa) == sign(fb)

error('Function must change sign on the interval')

end

c = a;

fc = fa;

d = b - c;

e = d;

% Main loop, exit from middle of the loop

while fb ~= 0

% The three current points, a, b, and c, satisfy:

% f(x) changes sign between a and b.

% abs(f(b)) <= abs(f(a)).

% c = previous b, so c might = a.

% The next point is chosen from

% Bisection point, (a+b)/2.

% Secant point determined by b and c.

% Inverse quadratic interpolation point determined

% by a, b, and c if they are distinct.

if sign(fa) == sign(fb)

a = c; fa = fc;

d = b - c; e = d;

end

if abs(fa) < abs(fb)

c = b; b = a; a = c;

fc = fb; fb = fa; fa = fc;

end

% Convergence test and possible exit

m = 0.5*(a - b);

tol = 2.0*eps*max(abs(b),1.0);

if (abs(m) <= tol) | (fb == 0.0)

break

end

% Choose bisection or interpolation

if (abs(e) < tol) | (abs(fc) <= abs(fb))

% Bisection

d = m;

e = m;

else

% Interpolation

s = fb/fc;

if (a == c)

% Linear interpolation (secant)

p = 2.0*m*s;

q = 1.0 - s;

else

% Inverse quadratic interpolation

q = fc/fa;

r = fb/fa;

p = s*(2.0*m*q*(q - r) - (b - c)*(r - 1.0));

q = (q - 1.0)*(r - 1.0)*(s - 1.0);

end;

if p > 0, q = -q; else p = -p; end;

% Is interpolated point acceptable

if (2.0*p < 3.0*m*q - abs(tol*q)) & (p < abs(0.5*e*q))

e = d;

d = p/q;

else

d = m;

e = m;

end;

end

% Next point

c = b;

fc = fb;

if abs(d) > tol

b = b + d;

else

b = b - sign(b-a)*tol;

end

fb = feval(F,b,varargin{:});

end

HARMONIC

function ydot = harmonic(t,y)

ydot = [y(2); -y(1)];

HARMONIC DAM PED

function ydot = harmonic_damped(t,y)

alpha = .5;

ydot = [y(2); -y(1) - alpha *y(2)];

HORNER

% This function evaluates

% the polynomial of degree n-1 with coefficients a1,...,an

% at the point z by using Horner's algorithm

function p = horner(a,z)

n = length(a);

m = length(z);

p = a(1)*ones(m,1);

for k = 2:n

p = z.*p + a(k);

end

LUTX

function [L,U,p] = lutx(A)

%LUTX Triangular factorization, textbook version

% [L,U,p] = lutx(A) produces a unit lower triangular matrix L,

% an upper triangular matrix U, and a permutation vector p,

% so that L*U = A(p,:)

[n,n] = size(A);

p = (1:n)';

for k = 1:n-1

% Find index of largest element below diagonal in k-th column

[r,m] = max(abs(A(k:n,k)));

m = m+k-1;

% Skip elimination if column is zero

if (A(m,k) ~= 0)

% Swap pivot row

if (m ~= k)

A([k m],:) = A([m k],:);

p([k m]) = p([m k]);

end

% Compute multipliers

i = k+1:n;

A(i,k) = A(i,k)/A(k,k);

% Update the remainder of the matrix

j = k+1:n;

A(i,j) = A(i,j) - A(i,k)*A(k,j);

end

% Separate result

L = tril(A,-1) + eye(n,n);

U = triu(A);

MGS

% Modified Gram-Schmidt algorithm (stable)

function [Q,R] = mgs(A)

[m,n] = size(A);

% Initializating variables

R = zeros(m,n);

S = zeros(m,m);

Q = zeros(m,m);

V = [A floor(10*rand(m,m-n))];

for i = 1:n

R(i,i)=norm(V(:,i));

Q(:,i)=V(:,i)/R(i,i);

for j = i+1:n

R(i,j) = Q(:,i)'*V(:,j);

V(:,j)= V(:,j)-R(i,j)*Q(:,i);

end

for i = n+1:m

S(i,i)=norm(V(:,i));

Q(:,i)=V(:,i)/S(i,i);

for j = i+1:n

S(i,j) = Q(:,i)'*V(:,j);

V(:,j)= V(:,j)-S(i,j)*Q(:,i);

end

ODE23TX

function [tout,yout] = ode23tx(F,tspan,y0,arg4,varargin)

%ODE23TX Solve non-stiff differential equations. Textbook version of ODE23.

% ODE23TX(F,TSPAN,Y0) with TSPAN = [T0 TFINAL] integrates the system

% of differential equations y' = f(t,y) from t = T0 to t = TFINAL. The

% initial condition is y(T0) = Y0. The input argument F is the

% name of an M-file, or an inline function, or simply a character string,

% defining f(t,y). This function must have two input arguments, t and y,

% and must return a column vector of the derivatives, y'.

% With two output arguments, [T,Y] = ODE23TX(...) returns a column

% vector T and an array Y where Y(:,k) is the solution at T(k).

% With no output arguments, ODE23TX plots the emerging solution.

% ODE23TX(F,TSPAN,Y0,RTOL) uses the relative error tolerance RTOL

% instead of the default 1.e-3.

% ODE23TX(F,TSPAN,Y0,OPTS) where OPTS = ODESET('reltol',RTOL, ...

% 'abstol',ATOL,'outputfcn',@PLOTFUN) uses relative error RTOL instead

% of 1.e-3, absolute error ATOL instead of 1.e-6, and calls PLOTFUN

% instead of ODEPLOT after each successful step.

% More than four input arguments, ODE23TX(F,TSPAN,Y0,RTOL,P1,P2,...),

% are passed on to F, F(T,Y,P1,P2,...).

% ODE23TX uses the Runge-Kutta (2,3) method of Bogacki and Shampine (BS23).

% Example

% tspan = [0 2*pi];

% y0 = [1 0]';

% F = '[0 1; -1 0]*y';

% ode23tx(F,tspan,y0);

% See also ODE23.

% Initialize variables.

rtol = 1.e-3;

atol = 1.e-6;

plotfun = @odeplot;

if nargin >= 4 & isnumeric(arg4)

rtol = arg4;

elseif nargin >= 4 & isstruct(arg4)

if ~isempty(arg4.RelTol), rtol = arg4.RelTol; end

if ~isempty(arg4.AbsTol), atol = arg4.AbsTol; end

if ~isempty(arg4.OutputFcn), plotfun = arg4.OutputFcn; end

end

t0 = tspan(1);

tfinal = tspan(2);

tdir = sign(tfinal - t0);

plotit = (nargout == 0);

threshold = atol / rtol;

hmax = abs(0.1*(tfinal-t0));

t = t0;

y = y0(:);

% Make F callable by feval.

if ischar(F) & exist(F)~=2

F = inline(F,'t','y');

elseif isa(F,'sym')

F = inline(char(F),'t','y');

end

% Initialize output.

if plotit

feval(plotfun,tspan,y,'init');

else

tout = t;

yout = y.';

end

% Compute initial step size.

s1 = feval(F, t, y, varargin{:});

r = norm(s1./max(abs(y),threshold),inf) + realmin;

h = tdir*0.8*rtol^(1/3)/r;

% The main loop.

while t ~= tfinal

hmin = 16*eps*abs(t);

if abs(h) > hmax, h = tdir*hmax; end

if abs(h) < hmin, h = tdir*hmin; end

% Stretch the step if t is close to tfinal.

if 1.1*abs(h) >= abs(tfinal - t)

h = tfinal - t;

end

% Attempt a step.

s2 = feval(F, t+h/2, y+h/2*s1, varargin{:});

s3 = feval(F, t+3*h/4, y+3*h/4*s2, varargin{:});

tnew = t + h;

ynew = y + h*(2*s1 + 3*s2 + 4*s3)/9;

s4 = feval(F, tnew, ynew, varargin{:});

% Estimate the error.

e = h*(-5*s1 + 6*s2 + 8*s3 - 9*s4)/72;

err = norm(e./max(max(abs(y),abs(ynew)),threshold),inf) + realmin;

% Accept the solution if the estimated error is less than the tolerance.

if err <= rtol

t = tnew;

y = ynew;

if plotit

if feval(plotfun,t,y,'');

break

end

else

tout(end+1,1) = t;

yout(end+1,:) = y.';

end

s1 = s4; % Reuse final function value to start new step.

end

% Compute a new step size.

h = h*min(5,0.8*(rtol/err)^(1/3));

% Exit early if step size is too small.

if abs(h) <= hmin

warning('Step size %e too small at t = %e.\n',h,t);

t = tfinal;

end

if plotit

feval(plotfun,[],[],'done');

end

ORBIT EULER

function orbit_euler(F,tspan,h,y0)

clc

%t0 = 0;

%h = 0.001;

%y0 = [1;0];

t0 = tspan(1);

tf = tspan(2);

t1 = t0 + h;

u = t0:h:tf;

v = u;

k = 1;

hs = h/10;

while t1 <= tf

y = euler(F,t0,y0,hs,t1);

u(k) = y(1);

v(k) = y(2);

plot(u(k),v(k),'b')

hold on

t0 = t1;

t1 = t0 + h;

y0 = y;

k = k+1;

end

PCHIPSLOPES

function d = pchipslopes(h,delta)

% PCHIPSLOPES Slopes for shape-preserving Hermite cubic

% interpolation. pchipslopes(h,delta) computes d(k) = P'(x(k)).

% Slopes at interior points

% delta = diff(y)./diff(x).

% d(k) = 0 if delta(k-1) and delta(k) have opposites signs

% or either is zero.

% d(k) = weighted harmonic mean of delta(k-1) and delta(k)

% if they have the same sign.

n = length(h)+1;

d = zeros(size(h));

k = find(sign(delta(1:n-2)).*sign(delta(2:n-1)) > 0) + 1;

w1 = 2*h(k)+h(k-1);

w2 = h(k)+2*h(k-1);

d(k) = (w1+w2)./(w1./delta(k-1) + w2./delta(k));

% Slopes at endpoints

d(1) = pchipendpoint(h(1),h(2),delta(1),delta(2));

d(n) = pchipendpoint(h(n-1),h(n-2),delta(n-1),delta(n-2));

% -------------------------------------------------------

PCHIPTX

function v = pchiptx(x,y,u)

%PCHIPTX Textbook piecewise cubic Hermite interpolation.

% v = pchiptx(x,y,u) finds the shape-preserving piecewise cubic

% interpolant P(x), with P(x(j)) = y(j), and returns v(k) = P(u(k)).

% See PCHIP, SPLINETX.

% First derivatives

h = diff(x);

delta = diff(y)./h;

d = pchipslopes(h,delta);

% Piecewise polynomial coefficients

n = length(x);

c = (3*delta - 2*d(1:n-1) - d(2:n))./h;

b = (d(1:n-1) - 2*delta + d(2:n))./h.^2;

% Find subinterval indices k so that x(k) <= u < x(k+1)

k = ones(size(u));

for j = 2:n-1

k(x(j) <= u) = j;

end

% Evaluate interpolant

s = u - x(k);

v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k)));

% -------------------------------------------------------

function d = pchipslopes(h,delta)

% PCHIPSLOPES Slopes for shape-preserving Hermite cubic

% interpolation. pchipslopes(h,delta) computes d(k) = P'(x(k)).

% Slopes at interior points

% delta = diff(y)./diff(x).

% d(k) = 0 if delta(k-1) and delta(k) have opposites signs

% or either is zero.

% d(k) = weighted harmonic mean of delta(k-1) and delta(k)

% if they have the same sign.

n = length(h)+1;

d = zeros(size(h));

k = find(sign(delta(1:n-2)).*sign(delta(2:n-1)) > 0) + 1;

w1 = 2*h(k)+h(k-1);

w2 = h(k)+2*h(k-1);

d(k) = (w1+w2)./(w1./delta(k-1) + w2./delta(k));

% Slopes at endpoints

d(1) = pchipendpoint(h(1),h(2),delta(1),delta(2));

d(n) = pchipendpoint(h(n-1),h(n-2),delta(n-1),delta(n-2));

% -------------------------------------------------------

function d = pchipendpoint(h1,h2,del1,del2)

% Noncentered, shape-preserving, three-point formula.

d = ((2*h1+h2)*del1 - h1*del2)/(h1+h2);

if sign(d) ~= sign(del1)

d = 0;

elseif (sign(del1) ~= sign(del2)) & (abs(d) > abs(3*del1))

d = 3*del1;

end

PESOSNC

% Este fichero almacena los

% pesos de las reglas de

% cuadratura de Newton-Cotes

% con m puntos

% m es un entero que satisface 2 <= m <= 11

% p es un vector columna que contiene

% los pesos para la regla de cuadratura

% de Newton-Cotes con m puntos.

function p = pesosNC(m);

if m==2

p=[1 1]'/2;

elseif m==3

p=[1 4 1]'/6;

elseif m==4

p=[1 3 3 1]'/8;

elseif m==5

p=[7 32 12 32 7]'/90;

elseif m==6

p=[19 75 50 50 75 19]'/288;

elseif m==7

p=[41 216 27 272 27 216 41]'/840;

elseif m==8

p=[751 3577 1323 2989 2989 1323 3577 751]'/17280;

elseif m==9

p=[989 5888 -928 10496 -4540 10496 -928 5888 989]'/28350;

elseif m==10

p=[2857 15741 1080 19344 5778 5778 19344 1080 15741 2857]'/89600;

else

p=[16067 106300 -48525 272400 -260550 427368 -260550 272400 -48525 106300 16067]'/598752;

end;

PIECELIN

function v = piecelin(x,y,u)

%PIECELIN Piecewise linear interpolation.

% v = piecelin(x,y,u) finds the piecewise linear L(x)

% with L(x(j)) = y(j) and returns v(k) = L(u(k)).

% First divided difference

delta = diff(y)./diff(x);

% Find subinterval indices k so that x(k) <= u < x(k+1)

n = length(x);

k = ones(size(u));

for j = 2:n-1

k(x(j) <= u) = j;

end

% Evaluate interpolant

s = u - x(k);

v = y(k) + s.*delta(k);

POLYINTERP

function v = polyinterp(x,y,u)

%POLYINTERP Polynomial interpolation.

% v = POLYINTERP(x,y,u) computes v(j) = P(u(j)) where P is the

% polynomial of degree d = length(x)-1 with P(x(i)) = y(i).

% Use Lagrangian representation.

% Evaluate at all elements of u simultaneously.

n = length(x);

v = zeros(size(u));

for k = 1:n

w = ones(size(u));

for j = [1:k-1 k+1:n]

w = (u-x(j))./(x(k)-x(j)).*w;

end

v = v + w*y(k);

end

QRSTEPS

function [A,b] = qrsteps(A,b)

%QRSTEPS Orthogonal-triangular decomposition.

% Demonstrates M-file version of built-in QR function.

% R = QRSTEPS(A) is the upper trapezoidal matrix R that

% results from the orthogonal transformation, R = Q'*A.

% With no output argument, QRSTEPS(A) shows the steps in

% the computation of R. Press <enter> key after each step.

% [R,bout] = QRSTEPS(A,b) also applies the transformation to b.

[m,n] = size(A);

if nargin < 2, b = zeros(m,0); end

% Householder reduction to triangular form.

if nargout == 0

clc

if ~isempty(b), b, end

pause

end

for k = 1:min(m-1,n)

% Introduce zeros below the diagonal in the k-th column.

% Use Householder transformation, I - rho*u*u'.

i = k:m;

u = A(i,k);

sigma = norm(u);

% Skip transformation if column is already zero.

if sigma ~= 0

if u(1) ~= 0, sigma = sign(u(1))*sigma; end

u(1) = u(1) + sigma;

rho = 1/(conj(sigma)*u(1));

% Update the k-th column.

A(i,k) = 0;

A(k,k) = -sigma;

% Apply the transformation to remaining columns of A.

j = k+1:n;

v = rho*(u'*A(i,j));

A(i,j) = A(i,j) - u*v;

% Apply the transformation to b.

if ~isempty(b)

tau = rho*(u'*b(i));

b(i) = b(i) - tau*u;

end

if nargout == 0

clc

if ~isempty(b), b, end

pause

end

QUADTX

function [Q,fcount] = quadtx(F,a,b,tol,varargin)

%QUADTX Evaluate definite integral numerically.

% Q = QUADTX(F,A,B) approximates the integral of F(x) from A to B

% to within a tolerance of 1.e-6. F is a string defining a function

% of a single variable, an inline function, a function handle, or a

% symbolic expression involving a single variable.

% Q = QUADTX(F,A,B,tol) uses the given tolerance instead of 1.e-6.

% Arguments beyond the first four, Q = QUADTX(F,a,b,tol,p1,p2,...),

% are passed on to the integrand, F(x,p1,p2,..).

% [Q,fcount] = QUADTX(F,...) also counts the number of evaluations

% of F(x).

% See also QUAD, QUADL, DBLQUAD, QUADGUI.

% Make F callable by feval.

if ischar(F) & exist(F)~=2

F = inline(F);

elseif isa(F,'sym')

F = inline(char(F));

end

% Default tolerance

if nargin < 4 | isempty(tol)

tol = 1.e-6;

end

% Initialization

c = (a + b)/2;

fa = feval(F,a,varargin{:});

fc = feval(F,c,varargin{:});

fb = feval(F,b,varargin{:});

% Recursive call

[Q,k] = quadtxstep(F, a, b, tol, fa, fc, fb, varargin{:});

fcount = k + 3;

% ---------------------------------------------------------

function [Q,fcount] = quadtxstep(F,a,b,tol,fa,fc,fb,varargin)

% Recursive subfunction used by quadtx.

h = b - a;

c = (a + b)/2;

fd = feval(F,(a+c)/2,varargin{:});

fe = feval(F,(c+b)/2,varargin{:});

Q1 = h/6 * (fa + 4*fc + fb);

Q2 = h/12 * (fa + 4*fd + 2*fc + 4*fe + fb);

if abs(Q2 - Q1) <= tol

Q = Q2 + (Q2 - Q1)/15;

fcount = 2;

else

[Qa,ka] = quadtxstep(F, a, c, tol, fa, fd, fc, varargin{:});

[Qb,kb] = quadtxstep(F, c, b, tol, fc, fe, fb, varargin{:});

Q = Qa + Qb;

fcount = ka + kb + 2;

end

RESTRICTED

function dydt = restricted(t,y)

mu = 1 / 82.45;

mustar = 1 - mu;

r13 = ((y(1) + mu)^2 + y(2)^2) ^ 1.5;

r23 = ((y(1) - mustar)^2 + y(2)^2) ^ 1.5;

dydt = [ y(3)

y(4)

(2*y(4) + y(1) - mustar*((y(1)+mu)/r13) - mu*((y(1)-mustar)/r23))

(-2*y(3) + y(2) - mustar*(y(2)/r13) - mu*(y(2)/r23)) ];

SIMPSON

% This function computes the

% Simpson quadrature rule for

% computing the integral of f(x) on

% the interval [a,b]

function I = simpson(fname,a,b)

fa = feval(fname,a);

fm = feval(fname, (a+b)/2);

fb = feval(fname, b);

I = ( fa + 4*fm + fb)*(b-a)/6;

SIMPSONCOM

% This function computes

% the Simpson composite quadrature rule

% with n subintervals

function I = simpsoncom(fname,a,b,n);

h = (b-a)/n;

x = [a:h:b];

f = feval(fname, x);

I = 0;

for k = 1:n

I = I + simpson('fname',x(k),x(k+1));

end

SPLINESLOPES

function d = splineslopes(h,delta)

% SPLINESLOPES Slopes for cubic spline interpolation.

% splineslopes(h,delta) computes d(k) = S'(x(k)).

% Uses not-a-knot end conditions.

% Diagonals of tridiagonal system

n = length(h)+1;

a = zeros(size(h)); b = a; c = a; r = a;

a(1:n-2) = h(2:n-1);

a(n-1) = h(n-2)+h(n-1);

b(1) = h(2);

b(2:n-1) = 2*(h(2:n-1)+h(1:n-2));

b(n) = h(n-2);

c(1) = h(1)+h(2);

c(2:n-1) = h(1:n-2);

% Right-hand side

r(1) = ((h(1)+2*c(1))*h(2)*delta(1)+h(1)^2*delta(2))/c(1);

r(2:n-1) = 3*(h(2:n-1).*delta(1:n-2)+h(1:n-2).*delta(2:n-1));

r(n) = (h(n-1)^2*delta(n-2)+(2*a(n-1)+h(n-1))*h(n-2)*delta(n-1))/a(n-1);

% Solve tridiagonal linear system

d = tridisolve(a,b,c,r);

SPLINTEX

function v = splinetx(x,y,u)

%SPLINETX Textbook spline function.

% v = splinetx(x,y,u) finds the piecewise cubic interpolatory

% spline S(x), with S(x(j)) = y(j), and returns v(k) = S(u(k)).

% See SPLINE, PCHIPTX.

% First derivatives

h = diff(x);

delta = diff(y)./h;

d = splineslopes(h,delta);

% Piecewise polynomial coefficients

n = length(x);

c = (3*delta - 2*d(1:n-1) - d(2:n))./h;

b = (d(1:n-1) - 2*delta + d(2:n))./h.^2;

% Find subinterval indices k so that x(k) <= u < x(k+1)

k = ones(size(u));

for j = 2:n-1

k(x(j) <= u) = j;

end

% Evaluate spline

s = u - x(k);

v = y(k) + s.*(d(k) + s.*(c(k) + s.*b(k)));

% -------------------------------------------------------

TRAPCOM

% This function computes

% the composite trapezoid quadrature rule with

% n subintervals

function I = trapezoid(fname,a,b,n);

h = (b-a)/n;

x = [a:h:b];

f = feval(fname, x);

I = h*(sum(f)-(f(1)+ f(length(f)))/2);

TWOBODY

function ydot = twobody(t,y)

r = sqrt(y(1)^2 + y(2)^2);

ydot = [y(3); y(4); -y(1)/r^3 ; -y(2)/r^3];

Tal vez te pueda interesar:

Descargar

Enviado por:	El remitente no desea revelar su nombre
Idioma:	castellano
País:	España

Palabras clave:

Cálculo matricial Matemáticas y computación Comandos Básicos Matemáticos Matrices Vectores Polinomios Derivadas y primitivas Teorema de los ceros de Bolzano Método de Newton

Te va a interesar