Usuario:Grupo 4

Trabajo realizado por estudiantes
Título	Fluido por encima de una placa plana. Grupo 4
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2015-16
Autores	Ruben Martos López, Guillermo Megino León, Silviu Popa Alejandro Sistac Ara,
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

2 Resolución numérica del problema.

2.1 Resolución mediante Euler modificado.

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
for j=1:G
   f0=[0 0 k(j)];
   f(:,1)=f0;
for i=1:N
      K1=w(t(i),f(:,i));
      t(i)=t(i)+h/2;
      z(:,i)=f(:,i)+K1*h/2;
      f(:,i+1)=f(:,i)+h*w(t(i),z(:,i));
end
   Y(j)=f(2,end);
end
plot(k,Y)
for q=1:G
   a(q)=abs(Y(q)-1);
end
l=min(a);
for q=1:G
   b(q)=abs(Y(q)-1);
if min(b)==l
      k(q)
      break
end

2.2 Resolución mediante Runge-Kutta de orden 4.

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
for j=1:G
   f0=[0 0 k(j)];
   f(:,1)=f0;
for i=1:N
      K1=w(t(i),f(:,i));
      t(i)=t(i)+h/2;
      z(:,i)=f(:,i)+K1*h/2;
      K2=w(t(i),z(:,i));
      v(:,i)=f(:,i)+K2*h/2;
      K3=w(t(i),v(:,i));
      c(:,i)=f(:,i)+K3*h;
      K4=w(t(i)-h/2,c(:,i));
      f(:,i+1)=f(:,i)+(h/6)*(K1+2*K2+2*K3+K4);
end
   Y(j)=f(2,end);
end
plot(k,Y)
for q=1:G
    a(q)=abs(Y(q)-1);
end
l=min(a);
for q=1:G
   b(q)=abs(Y(q)-1);
if min(b)==l
      k(q)
      break
end

2.3 Resolución mediante Euler.

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
for j=1:G
   f0=[0 0 k(j)];
   f(:,1)=f0;
for i=1:N
      f(:,i+1)=f(:,i)+h*w(t(i),f(:,i));
end
   Y(j)=f(2,end);
end
plot(k,Y)
for q=1:G
   a(q)=abs(Y(q)-1);
end
l=min(a);
for q=1:G
   b(q)=abs(Y(q)-1);
if min(b)==l
   k(q)
   q
   break
end
end

3 Interpretación de la derivada de la función de forma f(η)

clc
clear all
%Datos iniciales
%Numero de ecuaciones
ne=3;
%Tiempo inicial
t0=0;
%Tiempo final
tn=20;
%Tamaño del salto
h=0.05;
%Numero de elementos
N=round((tn-t0)/h);
%Determinacion del vector t
t=linspace(t0,tn,N+1);
%Vector k, valor inicial de f(3)
k0=0.1;
kf=1;
k=k0:0.01:kf;
G=length(k);
%Sistema de ecuaciones
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
%Resolucion
   f0=[0 0 k(24)];
   f(:,1)=f0;
   for i=1:N
  K1=w(t(i),f(:,i));
  t(i)=t(i)+h/2;
  z(:,i)=f(:,i)+K1*h/2;
  K2=w(t(i),z(:,i));
  v(:,i)=f(:,i)+K2*h/2;
  K3=w(t(i),v(:,i));
  c(:,i)=f(:,i)+K3*h;
  K4=w(t(i)-h/2,c(:,i));
  f(:,i+1)=f(:,i)+(h/6)*(K1+2*K2+2*K3+K4);

   end
 O=ones(1,N+1);
    plot(t,f(2,:))
 for j=1:N+1   
   if abs(f(2,j)-O)<0.01
       t(j)
       j
       break
   end
end

4 Estudio de la velocidad en función del factor de forma f(η)

4.1 Resolución numérica y representación gráfica de la componente U₁ del campo de velocidades

clc
clear all
 %Datos iniciales
ne=3;
h=0.05;
x=[0.05 0.2 0.4 0.6 0.8];
y0=0;
yn=3;
y=y0:h:yn;
P=length(y);
N=round((yn-y0)/h);
w=inline('[f(2);f(3);-1/2*f(1)*f(3)]','t','f');
f=zeros(ne,N+1);
for j=1:length(x)
   H2=h*sqrt(2/x(j));
   f0=[0 0 0.33];
   f(:,1)=f0;
   for i=1:N
   t(i)=y(i)*sqrt(2/x(j))+H2/2;
   K1=w(t(i),f(:,i));
   q(:,i)=f(:,i)+K1*H2/2;
   f(:,i+1)=f(:,i)+H2*w(t(i),q(:,i));
   end
   Z(j,:)=f(2,:);
   hold on 
   plot(y,2*Z(j,:))
   hold off
end

5 Representación gráfica e interpretación del borde de la capa limite

clc
clear all
x0=0;
xn=10;
h=0.05;
x=x0:h:xn;
eta=5.225;
for i=1:length(x)
    y(i)=eta*sqrt(x(i)/2);
end
plot(x,y);