Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

MTODOS DE UN PASO
EULER GENRICO
function [w]=euler(a, b, N, alfa) // a<x<b, N numero de subintervalos. deff('[q]=f(t,w)','q=forma de f') //definimos la forma de f segun corresp h=(b-a)/N t(1)=a w(1)=alfa

for i=1:N w(i+1)=w(i)+h*f(t(i),w(i)) t(i+1)=t(i)+h; end plot (t,w) endfunction

HEUN GENRICO
function [w] = heun(a,b,N,alfa) deff('[p]=f(t,w)','p=forma de f')//defino f h=(b-a)/N; w(1)=alfa t(1)=a

for i=1:N t(i+1)=t(i)+h; //paso temporal p(i)=f(t(i),w(i)) //defino p(i) como la pendiente en f(t(i),w(i)) w0(i+1)=w(i)+h*f(t(i),w(i)); wp(i+1)=f(t(i+1),w0(i+1)); //defino wp como "omega prima" w(i+1)=w(i)+(h/2)*(p(i)+wp(i+1)) end plot(t,w) endfunction

Ernesto E. Merlo, ING-237

Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

PUNTO MEDIO GENRICO

function [w] = pmedio(a,b,N,alfa) deff('[f]=f(t,w)','f=formadef')//defino la funcion f interna en el programa h=(b-a)/N; w(1)=alfa t(1)=a

for i=1:N t(i+1)=t(i)+h w(i+1)=w(i)+h*f(t(i)+h/2,w(i)+(h/2)*f(t(i),w(i))) end plot(t,w) endfunction

RUNGE-KUTTA ORDEN 4 GENRICO function [w]=rk4(a,b,N,alfa) deff('[f]=f(t,w)','f=formadef') //defino f para cada problema w(1)=alfa t(1)=a h=(b-a)/N //se debe abrir f como funcion adjunta for i=1:N

k1=h*f(t(i),w(i)) k2=h*f(t(i)+h/2,w(i)+k1/2) k3=h*f(t(i)+h/2,w(i)+k2/2) k4=h*f(t(i)+h,w(i)+k3) w(i+1)=w(i)+(1/6)*(k1+2*k2+2*k3+k4) t(i+1)=t(i)+h end plot(t,w) endfunction

Ernesto E. Merlo, ING-237

Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

EULER GENRICO GENERALIZADO PARA SISTEMAS (de 2 ecuaciones diferenciales lineales) function [w1,w2] = eulersistema(a,b,N,alfa1,alfa2) deff('[p]=f1(t,w1,w2)','p=formadef1')//definir f1 y f2 como corresponda segun el problema deff('[q]=f2(t,w1,w2)','q=formadef2') h=(b-a)/N t(1)=a //t varia en j w1(1)=alfa1 //w1 se corresponde con la aprox. a u1 variando en tiempo t(j) w2(1)=alfa2 //w2 se corresponde con la aprox. a u2 variando en tiempo t(j)

for j=1:N w1(j+1)=w1(j)+h*f1(t(j),w1(j),w2(j)) w2(j+1)=w2(j)+h*f2(t(j),w1(j),w2(j)) t(j+1)=t(j)+h end plot(t,w1) disp("se grafica w1 vs t por defecto.") endfunction

Ernesto E. Merlo, ING-237

Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

RUNGE-KUTTA ORDEN 4 GENRICO GENERALIZADO PARA SISTEMAS (lineal de 2 ED) function [w1,w2]=rksistema(a,b,N,alfa1,alfa2) deff('[p]=f1(t,w1,w2)','p=formadef1') //ser necesario definir cada f1 deff('[q]=f2(t,w1,w2)','q=formadef2') //ser necesario definir cada f2 h=(b-a)/N t(1)=a //t varia en j w1(1)=alfa1 w2(1)=alfa2 //cada wi varia en j disp("se graficar w1 por defecto ya que es el resultado") disp("mas util para ciertos problemas de orden superior") for j=1:N

t(j+1)=t(j)+h k11=h*f1(t(j),w1(j),w2(j)) k12=h*f2(t(j),w1(j),w2(j))// Kli, seria k sub ele corresp a Wi k21=h*f1(t(j)+(h/2),w1(j)+k11/2,w2(j)+k12/2) k22=h*f2(t(j)+(h/2),w1(j)+k11/2,w2(j)+k12/2) k31=h*f1(t(j)+(h/2),w1(j)+k21/2,w2(j)+k22/2) k32=h*f2(t(j)+(h/2),w1(j)+k21/2,w2(j)+k22/2) k41=h*f1(t(j+1),w1(j)+k31,w2(j)+k32) k42=h*f2(t(j+1),w1(j)+k31,w2(j)+k32) w1(j+1)=w1(j)+(1/6)*(k11+2*k21+2*k31+k41); //calculo los sucesivos a w2(j+1)=w2(j)+(1/6)*(k12+2*k22+2*k32+k42); //partir de los valores recien obtenidos

end plot(t,w1) endfunction

Ernesto E. Merlo, ING-237

Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

MTODOS DE DIFERENCIAS FINITAS

PARABOLICO EXPLICITO GENRICO
function [wt]=mpe(alfacuad, k, m, ti, tf, l, a, b) //alfacuad es alfacuadrado,k paso en t, h paso en x //tf tiempo final, ti tiempo inicial, l longitud del intervalo x //a valor de frontera inicial, b valor de frontera final deff('[f]=f(x)','f=formadef) //defino la funcion de condicion inicial h=l/m, landa=alfacuad*k/h^2, n=round((tf-ti)/k)

if landa<=0.5 then //Si se asegura la estabilidad, v1=(1-2*landa)*ones(m-1,1) v2=landa*ones(m-2,1) A=diag(v1,0)+diag(v2,1)+diag(v2,-1)

for i=1:m+1 //esta pensado as para incluir a los valores de frontera x(i)=(i-1)*h//forma de los x incluyendo los valores de frontera for j=1:n t=j*k w(i,1)=f(x(i))//asigno a cada valor lo que corresponde w(1,j)=a w(m+1,j)=b //con estos le digo que reemplacen las cond de front. w(2:m,j+1)=A*w(2:m,j)//la notacion 2:m me indica "desde el elemento 2 hasta //el elemento m" /Resuelve AWj=Wj+1 end end wt=w else disp("El mtodo ser inestable, por la eleccin de k (paso temporal)") wt=["no puede hallarse wij"] end endfunction

Ernesto E. Merlo, ING-237

Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

PARABOLICO IMPLICITO GENRICO

function [wt]=mpi(alfacuad, k, m, ti, tf, l, a, b) deff('[f]=f(x)','f=formadef') //f es condicion inicial h=l/m, landa=alfacuad*k/h^2, n=round(tf-ti)/k //defino variables

v1=(1+2*landa)*ones(m-1,1) v2=-landa*ones(m-2,1) A=diag(v1,0)+diag(v2,-1)+diag(v2,1) //formo la matriz A como suma de 3 diagonales

//definimos la matriz w que tendra en sus columnas los sucesivos wj for i=1:m+1 //esto esta pensado para incluir los valores de frontera x(i)=(i-1)*h //valores de x for j=1:n w(i,j)=f(x(i)) w(1,j)=a //aplica condicion de frontera inicial (fila 1 corresponde a x1) w(m+1,j)=b //aplica condicion de frontera final (fila 2 corresponde a xm+1) w(2:m,j+1)=inv(A)*w(2:m,j) //la notacion 2:m me indica "desde la componente //2 hasta la componente m". /resuelvo el sist AWj=Wj-1 end end wt=w

endfunction

Ernesto E. Merlo, ING-237

Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

METODO DE CRACK-NICOLSON GENRICO

function [wt]=mcn(alfacuad, k, m, ti, tf, l, a, b) deff('[f]=f(x)',['f=formadef]) //definir la forma de f cada ejercicio h=l/m landa=k*alfacuad/h^2 n=(tf-ti)/k

v1=(1+landa)*ones(m-1,1) v3=(1-landa)*ones(m-1,1) v2=(landa/2)*ones(m-2,1) v4=(-landa/2)*ones(m-2,1) A=diag(v4,-1)+diag(v1,0)+diag(v4,1) //formo A B=diag(v2,-1)+diag(v3,0)+diag(v2,1) //formo B R=inv(A) for i=1:m+1 x(i)=(i-1)*h//esto esta pensado para definir las x incluyendo las fronteras for j=1:n w(i,1)=f(x(i)) w(1,j)=a w(m+1,j)=b//defino para todos los j w(m+1,1)=b//defino para el j=1 w(2:m,j+1)=R*B*w(2:m,j) //resuelve Wj+1=inv(A)*B*Wj end end wt=w

endfunction

Ernesto E. Merlo, ING-237

Pseudocdigos de programacin para SciLab para Mtodos Computacionales en Ingeniera II

METODO HIPERBOLICO GENRICO

function [wt]=mh(alfacuad, m, k, l, ti, tf, a, b) deff('[f]=f(x)','f=forma de f') deff('[g]=g(x)','g=forma de g') h=l/m, alfa=sqrt(alfacuad), landa=alfa*k/h, n=round((tf-ti)/k) //tener cuidado con los valores para probar, n>1 //formamos la matriz A de (m-1)x(m-1) elementos p=(1-landa^2) q=landa^2 v1=(2*(1-landa^2))*ones(m-1,1) v2=(landa^2)*ones(m-2,1) A=diag(v1,0)+diag(v2,1)+diag(v2,-1)

//ahora quiero formar las primeras dos columnas de Wij //que ser una matriz de (m+1,n) donde cada columna Wi en el tiempo j for i=2:m //este conteo define solo los internos, los de frontera se agregan despues x(i)=(i-1)*h w(i)=f(x(i)), w(1)=a, w(m+1)=b//con esto defino la primer columna de W (Wi1) end //con este end me aseguro de haber formado todos los w11

//ahora defino los wi2 for i=2:m //los externos valdran cero por defecto w(i,2)=(0.5*q*w(i-1))+(p*w(i))+(0.5*q*w(i+1))-(k*g(x(i))) //calcula wi2 ahora que tengo wi1 end for j=1:n w(2:m,j+2)=A*w(2:m,j+1)-w(2:m,j)//calculo los sucesivos variando j desde 1 end wt=w

endfunction

Ernesto E. Merlo, ING-237