function U=dirich(f1,f2,f3,f4,a,b,h,tol,max1) %Datos% - f1, f2, f3, f4 son las funciones en el contorno almacenadas como% cadenas de caracteres.% - a y b son los extremos superiores de los intervalos [0,a] y [0,b]% - h es el incremento% - tol es la tolerancia%Resultado% - U es la matriz, análoga a la de la Tabla 10.6, en la que se% almacena la solución numérica.%Inicialización de los parámetros y de U n=fix(a/h)+1;m=fix(b/h)+1;ave=(a*(feval(f1,0)+feval(f2,0))+b*(feval(f3,0)+feval(f4,0)))/(2*a+2*b);U=ave*ones(n,m); %Condiciones de contornoU(1,1:m)=feval(f3,0:h:(m-1)*h)';U(n,1:m)=feval(f4,0:h:(m-1)*h)';U(1:n,1)=feval(f1,0:h:(n-1)*h);U(1:n,m)=feval(f2,0:h:(n-1)*h);U(1,1)=(U(1,2)+U(2,1))/2;U(1,m)=(U(1,m-1)+U(2,m))/2;U(n,1)=(U(n-1,1)+U(n,2))/2;U(n,m)=(U(n-1,m)+U(n,m-1))/2; %Parámetro de sobrerrelajación w=4/(2+sqrt(4-(cos(pi/(n-1))+cos(pi/(m-1)))^2)); %Mejora de las aproximaciones err=1;cnt=0;while((err>tol)&(cnt<=max1)) err=0; for j=2:m-1 for i=2:n-1 relx = w*(U(i,j+1)+U(i,j-1)+U(i+1,j)+U(i-1,j)-4*U(i,j))/4; U(i,j)=U(i,j)+relx; if (err<=abs(relx)) err=abs(relx); end end end cnt=cnt+1;endU=flipud(U'); %U=dirich('f1','f2','f3','f4',4,4,0.5,0.01,10)[x t]=meshgrid(0:0.05:1.5,0:0.05:1.5);surf(x,t,U)

function U=finedif(f,g,a,b,c,n,m)

% Construcción de la solución numérica de utt(x,t)=c2uxx(x,t) en:% R={(x,t):0<=x<=a,0<=t<=b} con u(0,t)=0, u(a,t)=0 para:% 0<=t<=b y u(x,0)=f(x), ut(x,0)=g(x) para 0<=x<=a% Datos % f=u(x,0) dada como una cadena caracteres 'f'.% g=ut(x,o) dada como una cadena de caracteres 'g'.% a y b son los extremos superiores de los intervalos [0 a] y [0 b];% c es la constante de la ecuación de ondas.% n y m es el número de nodos en [0 a] y [0 b] respectivamente. % Resultados% U - es la matriz en la que se almacena la solución numérica. % Inicialización de los parámetros y de U h=a/(n-1);k=b/(m-1);r=c*k/h;r2=r^2;r22=r^2/2;s1=1-r^2;s2=2-2*r^2;U=zeros(n,m); %Cálculo de las dos primeras filas for i=2:n-1 U(i,1)=feval(f,h*(i-1)); U(i,2)=s1*feval(f,h*(i-1))+k*feval(g,h*(i-1))... +r22*(feval(f,h*i)+feval(f,h*(i-2)));end % Cálculo de las demás filasfor j=3:m for i=2:(n-1) U(i,j)=s2*U(i,j-1)+r2*(U(i-1,j-1)+U(i+1,j-1))-U(i,j-2); endendU=U';[x t]=meshgrid(0:0.1:1,0:0.1:1);z=U;meshc(x,t,z)

function U=forwdif(f,c1,c2,a,b,c,n,m) %Datos

% - f=u(x,0) almacenada como una cadena de caracteres 'f'% - c1=u(0,t) y c2=u(a,t)% - a y b son los extremos derechos de [0 a] y [0 b]% - c es la constante de la ecuación del calor% - n y m son el número de nodos en [0,a] y [0,b]%Resultados% - U es la matriz de aproximaciones. % Inicialización de los parámetros y de Uh=a/(n-1);k=b/(m-1);r=c^2*k/h^2;s=1-2*r;U=zeros(n,m); %Condiciones de contornoU(1,1:m)=c1;U(n,1:m)=c2; %Construcción de la primera fila de UU(2:n-1,1)=feval(f,h:h:(n-2)*h)'; %Construcción de las demás filas de U for j=2:m for i=2:n-1 U(i,j)=s*U(i,j-1)+r*(U(i-1,j-1)+U(i+1,j-1)); endendU=U';[x t]=meshgrid(0:0.1:1,0:0.01:0.1);z=U,mesh(x,t,z)