Trabajo Final De Matlab-diego Mendoza.docx

 METODO DE EULER --FUNCION function dydx=f(x,y) dydx=y.*(x*x-1);

--CODIGO function euler(f,a,b,n,y0) %datos %f=el nombre de la función como string %a

=limite inferior

%b

=limite superior%h =longitud del segmento

%y0

=f(a)

%x

=es el vector x

%n

=numero de segmentos

%resultados %y

=es el vector f(x)

h=(b-a)/n; n=n+1; y=zeros(n,1); x=zeros(n,1); x(1)=a; y(1)=y0; fprintf('

x

y

\n');

fprintf('===================\n'); fprintf('%10.6f%10.6f\n',x(1),y(1)); for i=1:n-1 x(i+1)=a+h*i; y(i+1)=y(i)+h*feval('f',x(i),y(i)); fprintf('%10.6f%10.6f\n',x(i+1),y(i+1)); end

//RESULTADO >> euler('x^2',0,1,10,0.1) x

y

=================== 0.000000 0.100000 0.100000 0.090000 0.200000 0.081090 0.300000 0.073305 0.400000 0.066635 0.500000 0.061037 0.600000 0.056459 0.700000 0.052846 0.800000 0.050151 0.900000 0.048345 1.000000 0.047427

 METODO DE RK4 --CODIGO function RK4(f,a,b,n,y0) %datos %f

=el nombre de la función como string

%a

=limite inferior

%b

=limite superior

%h

=longitud del segmento

%y0

=f(a)

%x

=es el vector x

%n

=numero de segmentos

%resultados %y

=es el vector f(x) h=(b-a)/n; n=n+1;

y=zeros(n,1);

x=zeros(n,1); x(1)=a; y(1)=y0; fprintf('

x

y

fprintf('====================\n'); fprintf('%10.6f%10.6f\n',x(1),y(1)); for i=1:n-1 k1=feval('f',x(i),y(i)); k2=feval('f',x(i)+h/2,y(i)+h*k1/2); k3=feval('f',x(i)+h/2,y(i)+h*k2/2); k4=feval('f',x(i)+h,y(i)+h*k3); x(i+1)=a+h*i; y(i+1)=y(i)+h*(k1+2*k2+2*k3+k4)/6; fprintf('%10.6f%10.6f\n',x(i+1),y(i+1)); end

--RESULTADO >> RK4('4*exp(0.8*x)-0.5*y',0,2,4,0.2) x

y

==================== 0.000000 0.200000 0.500000 0.126468 1.000000 0.102669 1.500000 0.137399 2.000000 0.386643

\n');

 METODO DE HEUN --CODIGO function heun(f,a,b,n,y0) %datos %f

=el nombre de la función como string

%a

=limite inferior

%b

=limite superior

%h

=longitud del segmento

%y0

=f(a)

%x

=es el vector x

%n

=numero de segmentos

%resultados %y

=es el vector f(x)

h=(b-a)/n; n=n+1; y=zeros(n,1); x=zeros(n,1); x(1)=a; y(1)=y0; p(1)=0; fprintf('

x

p

y

\n');

fprintf('===============================\n'); fprintf('%10.6f%10.6f%10.6f\n',x(1),p(1),y(1)); for i=1:n-1 p(i+1)=y(i)+h*feval('f',x(i),y(i)); x(i+1)=a+h*i; y(i+1)=y(i)+0.5*h*(feval('f',x(i),y(i))+feval('f',x(i+1),p(i+1))); fprintf('%10.6f%10.6f%10.6f\n',x(i+1),p(i+1),y(i+1)); end

--RESULTADO

>> heun('x^2',0,1,10,0.1) x

p

y

=============================== 0.000000 0.000000 0.100000 0.100000 0.090000 0.090545 0.200000 0.081581 0.082147 0.300000 0.074261 0.074825 0.400000 0.068016 0.068564 0.500000 0.062805 0.063329 0.600000 0.058579 0.059080 0.700000 0.055299 0.055779 0.800000 0.052934 0.053404 0.900000 0.051481 0.051954 1.000000 0.050966 0.051460

 APLICACION 01 --FUNCION 01 % p.m function f=p(t) f=1; end

--FUNCION 02 % q.m function f=q(t) f=0; end

--FUNCION 03

% r.m function f=r(t) f=t*(t-1); end

--CODIGO 01 function [T,X] = findiff(p,q,r,a,b,alpha,beta,n) % Solucion del problema de valor de frontera usando diferencias finitas % x(t)‟‟ = p(t)x‟(t)+q(t)x(t)+r(t) % x(a) = alpha, x(b) = beta % Entradas % p, q, r Nombres de las funciones % a,b Limites del intervalo [a,b] % alpha Valor frontera izquierdo …. beta Valor frontera derecho % n numero de pasos % Salida % T Vector de abscisas X Vector de ordenadas T = zeros(1,n+1); X = zeros(1,n-1); Va = zeros(1,n-2); Vb = zeros(1,n-1); Vc = zeros(1,n-2); Vd = zeros(1,n-1); h = (b - a)/n; for j=1:n-1, Vt(j) = a + h*j; end for j=1:n-1, Vb(j) = -h^2*feval(r,Vt(j)); end Vb(1) = Vb(1) + (1 + h/2*feval(p,Vt(1)))*alpha; Vb(n-1) = Vb(n-1) + (1 - h/2*feval(p,Vt(n-1)))*beta; for j=1:n-1, Vd(j) = 2 + h^2*feval(q,Vt(j)); end

for j=1:n-2, Va(j) = -1 - h/2*feval(p,Vt(j+1)); end for j=1:n-2, Vc(j) = -1 + h/2*feval(p,Vt(j)); end X = trisys(Va,Vd,Vc,Vb); T = [a,Vt,b]; X = [alpha,X,beta];

--CODIGO 02 % trisys.m function X = trisys(A,D,C,B) %--------------------------------------------------------% Solucion del sistema tridiagonal % Entradas % A vector sub diagonal % D vector diagonal % C vector super diagonal % B vector del lado derecho % Salida % X vector solucion %--------------------------------------------------------n = length(B); for k = 2:n, mult = A(k-1)/D(k-1); D(k) = D(k) - mult*C(k-1); B(k) = B(k) - mult*B(k-1); end X(n) = B(n)/D(n); for k = (n-1):-1:1, X(k) = (B(k) - C(k)*X(k+1))/D(k);

end --RESULTADO >> [T, X] = findiff ('p', 'q', 'r', 0,1,0, 0,4) T= Columns 1 through 4 0 0.2500 0.5000 0.7500 Column 5 1.0000 X= Columns 1 through 4 0 0.0176 0.0269 0.0210 Column 5 0

 PROPUESTO 03  y=dsolve('Dy=(y-x-1)^2+2','y(0)=1','x'); >> simplify(y) ans = x + tan(x) + 1 >>ezplot('x + tan(x) + 1')

 PROPUESTO 04 >> y=dsolve('D2y-19.62*y=0','y(0)=0','Dy(0)=0.25','x'); >> simplify(y) ans = (5*218^(1/2)*sinh((3*218^(1/2)*x)/10))/1308 >> ezplot('(5*218^(1/2)*sinh((3*218^(1/2)*x)/10))/1308',[-1,1])