URGENTE!!!!¿sabes como resolver una INTEGRAL DOBLE? me das un ejemplo POR FAVOR!!!!!?

Answer 1

JAJAJA, LO QUE PEDIS ES MEDIO DENSO, PERO QUE SERIA DE YAHOO RESPUESTAS SI NO APARECE UN MAGO PARA DARTE LA RESPUESTA NO?

Te doy la solucion para una pc, ya que supongo que asi podes ver como funciona, aunque ya vi lo demas creo que vale un punto diferente, al menos aca esto es resolucion de una INTEGRAL DOBLE segun versa para una catedra, obviamente es mas facil como otros dicen, pero con esto tenes de todo para impresionar, fijate bien y luego decidis que abrevias o no.
Yo te paso la resolucion mas cientifica,.....vos decidis como la aplicas o no

ESPERO QUE CON ESTO ME LIGUE MIS 10 PUNTOS NO?


Este capítulo se basa en los ejemplos. Se propondrán una serie de ejercicios, se comentarán los algoritmos necesarios para resolverlos y finalmente se dará la solución a cada uno de ellos.


Metodología de programación
Quizás el mayor punto fuerte de Matlab o Octave no es su sencillez o su productividad sino la gran cantidad de funciones útiles que nos proporcionan. El núcleo de Matlab son sólo unos pocos Mbytes, en cambio la instalación completa ocupa tres CD-ROMs. Todo son funciones y archivos cuya razón de ser es hacernos la vida mucho más fácil.

Esto genera una ley un poco paradójica: cuanto menos programemos mejores serán nuestros scripts. Uno puede pensar que el código que escribe uno mismo es mejor, que si al intérprete le pasamos todas las líneas que debe usar irá más rápido. Que si tenemos que cargar más funciones en memoria haremos que la ejecución vaya más lenta... Es justo lo contrario. La mayoría de las funciones en Octave, por ejemplo, son funciones escritas en C++ o en Fortran compiladas de una manera especial para que Octave las interprete. Significa que la velocidad de ejecución de estas subrutinas será mucho más cercana a la velocidad de los programas binarios que a la del intérprete.


8.1 Ejercicio. Cálculo de un gradiente
Un cálculo muy común con una muestra de datos bidimensional es calcular una especie de gradiente numérico definido como:
grad(M)= 漂R> 缂R> 缂R> 缂R> 缂R> 缂R> 缂R> 輯FONT> M(i,j+1)-M(i,j-1)


2

M(i+1,j)-M(i-1,j)


2

?> ?> ?> ?> ?> ?> ?> ?NT> = 漂R> 輯FONT> δx′(M)
δy′(M)
?> ?NT>

Que genera una matriz de menor tamaño al perder los extremos. Escribir la función que devuelve el resultado en una variable tipo estructura: grad.x=δx′(M) y grad.y=δy′(M).


8.1.1 Guía para la resolución del ejercicio
Lo más complicado de este ejercicio es entender el enunciado. En el fondo es exactamente el mismo cálculo que un gradiente numérico con la diferencia de que en este caso no nos importa la distancia entre los puntos. Lo más práctico para calcular la fórmula es construir una submatriz por cada elemento de la ecuación en diferencias; para calcular la submatriz es necesario saber que el último elemento de cualquier matriz o vector se puede llamar end. Por ejemplo, si nosotros queremos los elementos del quinto al último del vector u pero no sabemos sus dimensiones, en vez de calcularlo lo más práctico es hacer > > u(5:end).


8.1.2 Solución del ejercicio
function[grad]=mygradient(M)

grad.x = 0.5*(M(:,3:end)-M(:,1:end-2));

grad.y = 0.5*(M(3:end,:)-M(1:end-2,:));
8.2 Ejercicio. Diseño de una tobera
Queremos diseñar una tobera convergente divergente. Para ello impondremos que el radio de salida sea 10 veces mayor que el radio de la garganta, y que la tobera se forma mediante dos parábolas, una con eje en y y otra con eje en x. Las condiciones serán que en el punto de empalme de los dos arcos haya continuidad tanto de la función como de la derivada primera. Las dos funciones a ajustar serán entonces (con el problema convenientemente adimensionalizado)

y-=Px2+0.1 y+=Ax+B
Entonces tenemos el sistema de ecuaciones siguiente, donde l es el punto en x donde se empalman los dos arcos:
2Pl= A


2Al+B


Pl2+0.1=Ax+B
A+B=1
Donde se demuestra que existe solución para P aproximadamente 0.9 Resolver el sistema de ecuaciones anterior y representar las soluciones de P=0.9, P=1 y P=1.2
8.2.1 Guía para la resolución del ejercicio
Lo primero que debemos hacer es entender el enunciado, para ver más claro lo que se pide, es bueno dibujar las dos curvas que se nos proponen, para ello haremos;
octave:1> fplot('[x. ^2+0.1,sqrt(x)]',[0,1])
Para resolver numéricamente sistemas de ecuaciones lineales usaremos la función fsolve. Para saber cómo funciona teclearemos en la línea de comandos:
octave:2> help fsolve
8.2.2 Solución del ejercicio
El script que nos resuelve el problema es el siguiente:
function[f]=ttb(x)
global Pi
f(1)=x(1)/(2*sqrt(x(1)*x(3)+x(2)))-2*Pi*x(3);
f(2)=Pi*x(3)**2+0.1-sqrt(x(1)*x(3)+x(2));
f(3)=sqrt(x(1)+x(2))-1;
end

function[f]=tobera(x,a,b,l,P)
if x f=P*(x*x)+0.1;
else
f=sqrt(a*x+b);
endif
end
x0=[1 1 1];
P=[0.9 1 1.2];
hold on
for i=1:3
global Pi
Pi=P(i);
[res]=fsolve('ttb',x0);
xcoord=linspace(0,1,100);
for j=1:100
tob(j)=tobera(xcoord(j),res(1),res(2),res(3),P(i));
end
plot(xcoord,tob)
end

hold off
El resultado lo vemos en la figura 8.1.

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



Figure 8.1: Resultado del script



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

8.3 Ejercicio. El atractor de Lorentz
Se quiere integrar la ecuación diferencial del Atractor de Lorentz, de ecuaciones:
x=a(y-x)
y=x(b-z)-y
z=xy-cz
con a=10, b=28, c= 8


3


y representarla en el espacio.


8.3.1 Guía para la resolución del Ejercicio
Primero escribimos la función correspondiente a la ecuación diferencial. En las subrutinas de resolución de EDOs siempre tenemos que introducir la función de la forma dx/dt=f(x,t) pudiendo ser cualquiera de las variables un vector.
La rutina de Octave que nos integra la ecuación diferencial es lsode; en Matlab utilizaremos ode45 porque el problema no es stiff. Utilizar la ayuda para saber de qué manera tenemos que escribir la función de la ecuación diferencial.
El comando para representar curvas paramétricas en 3 dimensiones es plot3.
Escribiremos todo en un script y lo ejecutaremos
8.3.2 Solución del ejercicio
8.3.2.1 Octave
x=0;
function xdot=func(x,t)
a=10;b=28;c=8/3;
xdot(1,1)=a*(x(2)-x(1));
xdot(2,1)=x(1)*(b-x(3))-x(2);
xdot(3,1)=x(1)*x(2)-c*x(3);
end

x0=[1;1;1];
t=linspace(0,50,5000);
tic;x=lsode( "func ",x0,t);toc
plot3(x(:,1),x(:,2),x(:,3))
dando como resultado la figura 8.2:

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



Figure 8.2: Resultado del script



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

y el tiempo necesario para resolver la EDO que es de 3.1488 segundos

La función lsode es mucho más versátil que las funciones para la integración de ecuaciones diferenciales de Matlab. También lo es en el sentido que es menos rigurosa en la manera de escribir la función, en este caso hemos escrito la variable xdot en forma de vector columna como nos obliga Matlab pero aceptaría igualmente el uso de un vector columna.

Octave permite definir la función en el script de modo que podemos resolver el problema con un único archivo.


8.3.2.2 Octave no stiff
El método de integración por defecto es explícito de modo que podemos acelerar el proceso si utilizamos un esquema Adams-Bashforth incluyendo este comando en el script:
lsode_options('integration method','non-stiff')
Podemos hacerlo sin miedo porque el problema del atractor de Lorentz, aunque es caótico, no es stiff. El tiempo de ejecución desciende a 2.0016 segundos.


8.3.2.3 Octave y C++
Ya hemos visto que podemos utilizar funciones escritas en otros lenguajes para aumentar la velocidad de nuestros scripts. En el caso de este ejercicio la velocidad conseguida por la función lsode es aceptable pero podemos rebajarla bastante más escribiendo la función de la ecuación diferencial en C++ y utilizando la aplicación mkoctfile. Ya hemos aprendido a hacerlo en la sección 7.4. Supongamos que ya disponemos del archivo eqlorentz.oct que usaremos del mismo modo que cualquier función. El script que resuelve el problema es:
t=linspace(0,50,5000);
tic;x=lsode( "eqlorentz ",[1;1;1],t);toc
plot3(x(:,1),x(:,2),x(:,3))
La nueva versión del script es capaz de resolver el problema en tan sólo 0.28444 segundos que es el 10% del tiempo anterior.


8.3.2.4 Octave y C++ no stiff
La máxima optimización se consigue de este modo. El tiempo se rebaja hasta 0.17314 segundos. Aunque todas estas consideraciones sobre velocidad pueden parecer inútiles debemos tener en cuenta que la integración de EDOs y EDPs son los problemas numéricos más exigentes en lo que respecta a uso de CPU. Entrar en este tipo de discusiones a menudo comporta grandes beneficios aunque en este caso sean irrisorios. Este último tiempo es comparable al obtenido con código enteramente escrito en C++ y sólo un poco más lento que el escrito en C o Fortran y el esfuerzo necesario para escribir el programa ha sido mucho menor.1


8.3.2.5 Octave, C++ y Fortran
Fortran es el lenguaje del cálculo matricial por excelencia. Hemos visto ya la manera de producir wrappers para funciones en Fortran. Su utilidad, más que para escribir pequeñas funciones donde el wrapper puede ser más largo que el código útil, es la de acoplar subrutinas que hagan uso de grandes cantidades de memoria con distintas precisiones en coma flotante. La velocidad de Fortran es aproximadamente la misma que la de C++ pero aporta ventajas distintas a la velocidad por lo dicho anteriormente.


8.3.2.6 Matlab
En Matlab necesitaremos dos archivos. El archivo de función devuelve el resultado en forma de vector columna (lorentz.m):
function xdot=lorentz(t,x)
a=10;b=28;c=8/3;
xdot(1,1)=a*(x(2)-x(1));
xdot(2,1)=x(1)*(b-x(3))-x(2);
xdot(3,1)=x(1)*x(2)-c*x(3);
end
Fijémonos que el orden de las variables de la cabecera x y t cambian según la rutina que usemos. El script será:
x=0;t=0;
tic;[t,x]=ode45(@lorentz,[0,50],[1,1,1]);toc
plot3(x(:,1),x(:,2),x(:,3))



8.4 Ejercicio. CÁLCULO DE UNA INTEGRAL DOBLE


La integral
I= ∫ ∞



-∞
∫ ∞



-∞
e-(x2+y2)dx dy

tiene como resultado π. Calcular la solución con un único comando.


8.4.1 Guía para la resolución del ejercicio
La función que implementa la integral doble en matlab es dblquad mientras que en Octave es quad2dg. Ambas funciones tienen una cabecera idéntica, las únicas diferencias son el nombre y el algoritmo de cálculo. Para pasarles la función a integrar como argumento usaremos una sentencia inline o con una función anónima. Las rutinas de integración tienen muchos problemas con las singularidades porque trabajan con una precisión limitada, esto nos impide usar límites de integración demasiado grandes y ni mucho menos podremos usar Inf . Con [-10,10]כ-10,10] vamos a conseguir 5 cifras significativas.

Por si alguien aún no está trabajando con Octave esta es la ayuda de la función quad2dg.


>> help quad2dg
quad2dg is the user-defined function from the file
/usr/share/octave/2.1.64/site/m/octave-forge/integration/quad2dg.m
usage: int = quad2dg('Fun',xlow,xhigh,ylow,yhigh)
or
int = quad2dg('Fun',xlow,xhigh,ylow,yhigh,tol)
This function is similar to QUAD or QUAD8 for 2-dimensional integration,
but it uses a Gaussian quadrature integration scheme.
int -{}- value of the integral
Fun -{}- Fun(x,y) (function to be integrated)
xlow -{}- lower x limit of integration
xhigh -{}- upper x limit of integration
ylow -{}- lower y limit of integration
yhigh -{}- upper y limit of integration
tol -{}- tolerance parameter (optional)

Note that if there are discontinuities the region of integration
should be broken up into separate pieces. And if there are singularities,
a more appropriate integration quadrature should be used
(such as the Gauss-Chebyshev for a specific type of singularity).
8.4.2 Solución del ejercicio
Tenemos múltiples opciones:
Con inline
Matlab:


>> dblquad(inline('exp(-(x.*x+y.*y))'),-10,10,-10,10)

ans = 3.1416

Octave:




>> quad2dg(inline('exp(-(x.*x+y.*y))'),-10,10,-10,10)

ans = 3.1416

Con una función anónima:
Matlab:




>> dblquad(@(x,y) exp(-(x. ^2+y. ^2)),-10,10,-10,10)

ans = 3.1416




Octave:




>> quad2dg(@(x,y) exp(-(x. ^2+y. ^2)),-10,10,-10,10)

ans = 3.1416



8.5 Ejercicio. Resolución de la ecuación de Laplace en un dominio bidimensional
Resolver la ecuación de Laplace en un dominio rectangular por diferencias finitas con condiciones de contorno Dirichlet. La ecuación de Laplace es la ecuación elíptica más sencilla: ∇2φ=0, que formulada en dos dimensiones se convierte en:
∂2φ


∂ x2
+ ∂2φ


∂ y2
=0

Si se toman diferencias finitas centradas de segundo orden para puntos dependiendo de i y j llegamos a la ecuación en diferencias siguiente:
φ(i+1,j)-2φ(i,j)+φ(i-1,)


Δ x2
+ φ(i,j+1)-2φ(i,j)+φ(i,j-1)


Δ y2
=0

Esta ecuación puede ser expresada en forma de sistema de ecuaciones lineales Aϕ=b, donde b es un término independiente que aparece al aplicar las condiciones de contorno y ϕ es la matriz de incógnitas φ expresada en forma de vector columna. La traslación del problema bidimensional a un vector de incógnitas es un paso que nos puede costar de entender pero si tenemos una matriz de incógnitas el tensor del sistema tendrá tres dimensiones con lo que ya no tendremos rutinas escritas para resolver el sistema.

Usaremos como parámetros dx=2, dy=3, n=50, m=50. No podemos utilizar muchos más elementos porque se nos quedaríamos sin memoria. Esta matriz de incógnitas se va a convertir en un vector nנm, lo que significa que la matriz del sistema va a tener (nנm)2 elementos. Para un dominio de 100 por 100 puntos llegamos a 108 puntos. Sería un buen caso de aplicación de matrices sparse


8.5.1 Guía para la resolución del ejercicio
Escribir una función place que implemente la transformación de la matriz al vector. La entrada serán los índices i y j y el número de elementos por columna n. La salida de la función será la posición en el vector posterior: place=i+n(j-1).
Crear la matriz del sistema con la ecuación en diferencias y la función creada en el apartado anterior
Poner las condiciones de contorno al gusto en el vector del término independiente y resolver el sistema lineal
Para resolver el sistema lineal del modo usual basta con hacer A.
Para resolver el sistema con matrices sparse primero creamos la matriz sparse con:
spA=sparse(A)
y luego resolveremos el sistema del modo usual. Es una buena idea eliminar la matriz del sistema de la memoria.
8.5.2 Solución del ejercicio
Primero definimos las constantes del problema
dx=2;dy=3;n=50;m=25;
A continuación creamos una función que reordena cualquier elemento de una matriz bidimensional en un vector en el que se concatenan las columnas.
function place=place(i,j,n)

place=i+n*(j-1);

end
Ahora definimos la matriz del sistema teniendo en cuenta que la incógnita no va a ser una matriz de dos dimensiones sino un vector. Esto hace que la matriz del sistema tenga en realidad nmנnm elementos.
A=zeros(n*m,n*m);
for i=1:n*m
A(i,i)=1;
end
for i=2:n-1
for j=2:m-1
A(place(i,j,n),place(i,j,n))=-2*(1/(dx*dx)+1/(dy*dy));
A(place(i,j,n),place(i-1,j,n))=1/(dx*dx);
A(place(i,j,n),place(i+1,j,n))=1/(dx*dx);
A(place(i,j,n),place(i,j-1,n))=1/(dy*dy);
A(place(i,j,n),place(i,j+1,n))=1/(dy*dy);
end
end
Una vez definida la matriz del sistema creamos el vector del término independiente que contiene las condiciones de contorno.
i=1;
for j=1:m
b(place(i,j,n))=sin(pi*(j-1)/(m-1));
end
i=n;
for j=1:m
b(place(i,j,n))=-sin(pi*(j-1)/(m-1));
end
j=1;
for i=1:n
b(place(i,j,n))=sin(pi*(i-1)/(n-1));
end
j=n;
for i=1:n
b(place(i,j,n))=0;
end
Y finalmente resolvemos el sistema lineal.
T=A \b';
T=reshape(T,n,m);
contour(T);
El resultado del script son las figuras 8.3 y 8.4.

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




Figure 8.3: Superficie solución



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

En el capítulo dedicado a la representación gráfica de soluciones hemos hablado sobre lo poco útiles que suelen ser las superfícies coloreadas para dar un resultado. Aunque una de estas representaciones pueda parecer muy explicativa e intuitiva nunca debemos perder de vista que lo que nos importa es mostrar un resultado y una superfície no lo consigue. Por muy bien que coloquemos los colores o las barras no conocemos el valor de cada punto con precisión y si encima colocamos líneas de nivel dotamos al gráfico de un exceso de información.

Aunque sea mucho más espartano y pueda parecer inapropiado la solución al problema son las curvas de nivel, en inglés contours . Aunque en un principio cueste un poco más entender la solución estamos ofreciendo muchísima más información de un modo netamente más simple. Se podría decir que no hay motivos para no utilizar las curvas de nivel a parte de intentar demostrar nuestra pericia en el uso de Matlab.

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




Figure 8.4: Superficie solución



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

Los gráficos de curvas de nivel pueden ajustarse perfectamente a nuestras necesidades. Una práctica muy común y beneficiosa es la de embeber el valor de cada curva dentro de la misma para no obligar al lector a mirar la leyenda.


8.5.2.1 Mejorar la solución. Cálculo de tiempos.
Si analizamos el código de la solución es evidente que estamos utilizando pocas funciones de la biblioteca. Además estamos haciendo una de las pocas cosas casi prohibidas en Matlab que es crear una matriz por fuerza bruta.

Un primer paso para mejorar la ejecución de un código es comprobar qué partes consumen más recursos e intentar mejorarlas aisladamente. Si no somos capaces de conseguirlo nos replantearemos la solución de un modo global. Para calcular los tiempos llenaremos el código con las funciones tic y toc.
dx=2;
dy=3;
n=50;
m=50;
function place=place(i,j,n)
place=i+n*(j-1);
end
A=zeros(n*m,n*m);
tic
for i=1:n*m
...
end
toc;disp('Creacion de la matriz del sistema'),disp(toc);tic;
i=1;
for j=1:m
b(place(i,j,n))=sin(pi*(j-1)/(m-1));
...

b(place(i,j,n))=0;
end
toc;disp('Creacion del vector b'),disp(toc);tic;
T=A \b';
toc;disp('Resolucion del sistema'),disp(toc);
T=reshape(T,n,m);
Los tiempos en cada paso son:
>> ejercicio4
Creacion de la matriz del sistema
1.0611
Creacion del vector b
0.017038
Resolucion del sistema
2.1457
Parece que debemos mejorar la construcción de la matriz del sistema y la resolución del mismo. ¿Cómo? El primer paso suele ser analizar la forma de la matriz con la función spy.
>> spy(A)
Que tiene como resultado el siguiente patrón (figura 8.5):

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




Figure 8.5: Patrón de elementos no nulos de la matriz del sistema



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

Esto nos demuestra que es una matriz n-diagonal por bandas lo que entra en la definición de matriz sparse. Parece que la solución a nuestros ligeros problemas de velocidad pueden solucionarse con el uso de este tipo de matrices.


8.5.2.2 Resolución del problema mediante matrices sparse(+)
Uno de los errores cometidos en la primera solución propuesta es que no utiliza las rutinas de creación de matrices disponibles. La alternativa es utilizar bucles for para crearla por fuerza bruta, solución de ningún modo recomendable. La mejor opción suele ser utilizar la función diag pero como estamos tratando con matrices sparse utilizaremos spdiags.


8.5.2.3 Resolución del problema con un método iterativo.
Para resolver el problema por un método directo es estrictamente necesario romper la matriz para convertirla en un vector. Esto hace que no podamos utilizar los operadores diferenciales en forma de matriz tal como veremos en el ejercicio 8.6. Un truco para mantener la integridad de la matriz y poder utilizar estos operadores es resolver el problema con un método iterativo. Aunque estemos multiplicando una matriz llena por matrices casi vacías el resultado se acelera considerablemente.

La teoría de este método es bastante sencilla. Si aplicamos el operador laplaciano a la matriz bidimensional el problema numérico se reescribe del siguiente modo.
Dxφ+φDy⊤=0

Este planteamiento carece de solución analítica pero sirve para plantear un problema iterativo en el que sólo hay que evaluar la expresión. Empezaremos la subrutina del mismo modo que la anterior, definiendo la función que recoloca los elementos de la matriz bidimensional.
dx=2;dy=3;global n=50;global m=50;
% resolucion del problema con un solver iterativo construccion de los
% operadores en diferencias
function place=place(i,j,n)
place=i+n*(j-1);
end
Ahora definiremos las matrices de los operadores diferenciales. Las definimos como variables globales al igual que las dimensiones de la matriz porque son necesarios dentro de la función que calcula el paso de iteración y ésta puede tener sólo un argumento de entrada.
global DX=diag([1,-2/(dx*dx).*ones(1,n-2),1],0).+...
diag([0,1/(dx*dx).*ones(1,n-2)],1).+...
diag([1/(dx*dx).*ones(1,n-2),0],-1);
global DY=diag([1,-2/(dy*dy).*ones(1,n-2),1],0).+...
diag([0,1/(dy*dy).*ones(1,n-2)],1).+...
diag([1/(dy*dy).*ones(1,n-2),0],-1);
A continuación la función que calcula el término en cada iteración. Notemos el hecho de que el argumento que recibe la función debe ser un vector por requerimientos de la rutina de resolución. Para calcular sobre las matrices utilizamos la función reshape.
function [AT]=calcstuff(T)
global n
global m
global DX
global DY
rT=reshape(T,n,m);
AT=reshape((DX*rT)+(rT*DY'),n*m,1);
end
Ahora es el turno del término independiente que se define gracias a la función place definida anteriormente.
i=1;j=1;b=1;
for j=1:m
b(place(i,j,n))=sin(pi*(j-1)/(m-1));
end
i=n;
for j=1:m
b(place(i,j,n))=-sin(pi*(j-1)/(m-1));
end
j=1;
for i=1:n
b(place(i,j,n))=sin(pi*(i-1)/(n-1));
end
j=n;
for i=1:n
b(place(i,j,n))=0;
end
Y finalmente resolvemos el sistema.
tic;T=pcg('calcstuff',b',tol=1e-6,maxit=100);
disp('Tiempo de calculo'),disp(toc)
El método utilizado para resolver el problema es el del ``Gradiente Conjugado Precondicionado''. Para utilizarlo debemos tener en cuenta que el número de iteraciones por defecto no es sufciente; lo corregimos y lo situamos en 100. El tiempo de proceso ha pasado de más de tres segundos a 0.249292. Este programa no sería posible sin el uso inteligente de las variables globales tal como las hemos utilizado en otros ejemplos.


8.6 Ejercicio. Un problema de calor unidimensional
Este ejercicio requiere la función trisolve que está disponible en Octave. Matlab no tiene ninguna rutina para resolver directamente sistemas tridiagonales porque trata este tipo de matrices como sparse. Para resolver el mismo problema en Matlab utilizaremos la función spdiags para montar la matriz sparse y luego resolveremos el sistema de ecuaciones del modo usual.

Este ejercicio es una pequeña demostración de hasta dónde puede llegar el ahorro de código con Matlab. Se trata de resolver la ecuación del calor unidimensional con condiciones de contorno Dirichlet. La discretización espacial serán diferencias finitas de sevundo órden y la temporal será un esquema Crank Nicholson. Tomaremos este esquema porque permite usar saltos de tiempo más grandes, de hecho nos permitirá llegar al resultado sin tener que calcular ninguna iteración temporal. La solución propuesta tiene sólo unas pocas líneas.

Una solución tan breve requiere algo de preparación matemática de modo que hay que estar muy atento al planteamiento analítico del problema.


Discretización de la ecuación completa
La ecuación del calor unidimensional es la EDP parabólica más sencilla:
∂tφ=∂xxφ
Como cualquier EDP parabólica se puede formular de la siguientes manera:
dφ


dt
=F(x,t)

Si utilizamos un esquema Crank Nicholson la discretización temporal será de la forma:
φn+1-φn


Δ t
= 1


2
( Fn+1+Fn )

Discretizando también el lado derecho de la ecuación con diferencias finitas de segundo orden llegamos a la ecuación discretizada:
φin+1- Δ t


2
漂R> 缂R> 缂R> 輯FONT> φi+1n+1-2φin+1+φi-1n+1


Δ x2
?> ?> ?> ?NT> =φin+ Δ t


2
漂R> 缂R> 缂R> 輯FONT> φn+1n-2φin+φi-1n


Δ x2
?> ?> ?> ?NT>



Formulación matricial del problema.
La ecuación anterior puede escribirse del modo siguiente:
φn+1- Δ t


2
Aφn+1=φn+ Δ t


2
Aφn

Agrupando términos, diciendo que B?=I?Δ t/2A y utilizando un Δ x=1 (en el fondo es sólo adimensionalizar):
B-φn+1=B+φn
donde la matriz A tiene la forma:
A= 漂R> 缂R> 缂R> 缂R> 缂R> 缂R> 缂R> 輯FONT> -2 1 0 0 ⋯ 0
1 -2 1 0 0
0 1 -2 1 0
⋮ ⋰ ⋮
0 0 1 -2 1
0 ⋯ 0 0 1 -2
?> ?> ?> ?> ?> ?> ?> ?NT>



Condiciones de contorno e iniciales.
Condiciones de contorno: φ=0 en x=0 y φ=1 en x=10
Condiciones iniciales: φ(x)=0 en los puntos interiores.
Dar una solución al problema para cualquier tiempo inferior a t=10.
8.6.1 Guía para la resolución del ejercicio
La dificultad de la resolución de este problema radica en la correcta situación de las condiciones de contorno en la matriz, el uso correcto de las submatrices y entender el funcionamiento de trisolve.


8.6.2 Solución del ejercicio (Octave)
Parte del ejercicio es entender qué estamos haciendo para llegar a la solución. Lo bueno de abreviar en la escritura es que además el código resultante suele ser mucho más rápido.
phi=zeros(10,1); phi(10,1)=1;
dt=0.1;
for i=1:50 #Activate loop in case of precision leak
phi(2:9)=(1-dt)*phi(2:9)+0.5*dt*phi(1:8)+0.5*dt*phi(3:10);
phi=trisolve(-[0.5*dt*ones(8,1);0],[1;(1+dt)*ones(8,1);1],...
-[0;0.5*dt*ones(8,1)],phi);
end
De este modo llegamos al resultado del problema con una precisión aceptable. La figura del perfil de temperaturas es:

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



Figure 8.6: Figura solución



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

8.6.3 Solución del ejercicio (Matlab)
Aunque la resolución de sistemas tridiagonales es un clásico del cálculo numérico Matlab no cuenta con ninguna función para ello. El motivo es que decidieron incluir las matrices tridiagonales dentro de las matrices sparse. Las matrices tridiagonales o con pocas bandas se resuelven eficientemente mediante un método directo con lo que en este caso Matlab tiene un error de diseño.

No tendría ningún sentido programarse una función que resolviera sistemas tridiagonales en Matlab porque todas estas rutinas son en el fondo interfaces a bibliotecas en Fortran y en C por motivos de velocidad.

En este caso la dimensión de la matriz no justifica el uso del almacenamiento sparse de modo que utilizaremos la función diag para crear una matriz según sus diagonales:
phi=zeros(10,1); phi(10,1)=1; dt=0.1;
B=diag([1,(1+dt)*ones(1,8),1],0)+...
diag(-[0.5*dt*ones(1,8),0],-1)+...
diag(-[0,0.5*dt*ones(1,8)],1);
for i=1:50
phi(2:9)=(1-dt)*phi(2:9)+0.5*dt*phi(1:8)+0.5*dt*phi(3:10);
phi=B \phi;
end
8.6.4 Comprobación de la evolución temporal
La siguiente modificación del programa permite comprobar la evolución de los pasos temporales sacando por pantalla varios estados temporales de la solución. Aunque el problema es estable para saltos temporales muy grandes preferiremos acortarlo por motivos de precisión. Representaremos la solución en los tiempos t=[0.1,1,2,5,10,20,50,100]. Muy atentos al modo de definir la condición lógica que indica cuándo pintar la curva:
phi=zeros(10,1); phi(10,1)=1;dt=0.1;
# times to output
tout=[0.1,1,2,5,10,20,50,100];
hold on
for i=1:500
phi(2:9)=(1-dt)*phi(2:9)+0.5*dt*phi(1:8)+0.5*dt*phi(3:10);
phi=trisolve(-[0.5*dt*ones(8,1);0],[1;(1+dt)*ones(8,1);1],...
-[0;0.5*dt*ones(8,1)],phi);
if sum((dt*i)*ones(1,length(tout))==tout)==1
plot(phi)
end
end
hold off
title('Solucion de la ecuacion del calor unidimensional')
xlabel('Longitud (adimensional)')
ylabel('Temperatura (adimensional)')
legend({'0.1','1','2','5','10','20','50','100'},2)
Llegamos finalmente a la siguiente gráfica:

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



Figure 8.7: Evolución del perfil de temperaturas



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


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

1
El uso de lenguajes ``pegamento'' es una práctica cada día más habitual. Escribir un código enteramente en C o Fortran es bastante costoso tanto por el tiempo de desarrollo como por el esfuerzo necesario. Muchas veces se empieza con un prototipo escrito con un lenguaje de RAD (Rapid Application Development) y posteriormente se reescriben sus partes.

El mayor obstáculo en los proyectos grandes es que deben participar varios desarrolladores. Cada uno puede tener sus preferencias y las guías de estilo no siempre son una solución efectiva porque sólo arreglan el formato del código escrito. Todas estas dificultades unidas al aumento de velocidad de lenguajes como Python o Matlab ha llevado a la creación del concepto de ``lenguaje pegamento''.

Se trata de dejar a medias el código entre el prototipo y el programa definitivo. Partiendo del prototipo se van analizando las partes que consumen más recursos y, sin cambiar las cabeceras ni el esquema de variables, se reescriben en algún lenguaje rápido. Se para cuando se consigue un equilibrio entre nivel de interactividad, velocidad y coste de mantenimiento de código.

Probablemente el mejor ``lenguaje pegamento'' sea Python gracias a Pyrex, SWIG y F2Py. El primero es capaz de convertir el código Python en C++ automáticamente (rendimiento máximo y coste cero) y el segudo y el tercero son generadores automáticos de interfaces, SWIG para C y C++ y F2Py para Fortran 77 y Fortran 95.
2
Este ejercicio se ha resuelto en un AMD Athlon 3000+.

Answer 2

Saludos!
El chiste es integrar primero lo de adentro y despues lo de mas afuera, por ejemplo si tomamos:
INT( INT( x+y-1))dxdy

integras primero con respecto a "x" tomando a "y" como constante y te queda:
=INT((x^2)/2 +xy-x)dy

ahora integras con respecto a "y" y te queda:
=(x^2)y/2+x(y^2)/2-xy

Espero que te halla servido este ejemplo

Answer 3

Si usted quiere tener éxito con el trading online, usted debe tener un historial probado o tiene que seguir un método preciso.
Si quieres empezar a ganar finalmente seria usted debe buscar en este sitio http://trading-on.info

Answer 4

jejeje suena dificl pero es simple... si te dicen doble integral es porq tienes 2 diferenciales (dx y dy) esto se utiliza para las áres .. si alguna ves te pregunta por integrla triple es porq tienes 3 diferenciales (dx,dy dz) se usa en volumnes... ahora ejemplo..

∫∫2dxdy = ∫2dx∫dy (así se ve menos fea)
despues inetgras:
(2∫dx)(∫dy) = (2x^2 / 2) (y)= (x^2) (y)
solución:
x^2y...Nota ("^"signifa que esta elevado a potencia)

No está dificil, depende mucho de los maestros...(no te desanimes si no le entiendes mucho al principio) suerte..

Answer 5

La mayor complicacion en resolver una integral doble no es la operacion en si (que se resuelve de forma similar a las simples, como te explicaron mas arriba), si no que es definir la region de integracion

acordate que para regiones mas complejas en las que una variable varia entre dos numeros, pero la otra varia entre dos funciones (en funcion de la otra variable) la integral que varia entre dos funciones se resuelve primero

Answer 6

Como en YR no se pueden escribir símbolo matemáticos, el signo de integral te lo escribiré con el corchete { y para señalarte los límite de la misma diré : < integral de a2 a a1 > por { de a2 a a1 siendo a2 el límite inferior y a1 el límite superior,,, Ejemplo. Hallar la integral doble definida. { de o a a { de ceo a (a^2-x^2 ^(1/2) de (x + y ) dy dx SOLUCIÓN:; Primero integras con respecto a y , considerando a dx como constante.{ de cero a “a” { de cero a raíz cuadrada de a cuadrada menos x cuadrada, multiplicado por (x + y ) dydx Como es una intergral directa queda:= { cero a “a” [ { de cero a raíz cuadra de a cuadrada menos x cuadrada (x+ y) dy] dx= { de cero a “a” [xy + y^2 /2 ] de cero a raiz cuadrada de a^2 – x^2 dx = { de cero a “a” (x por raíz cudrada de (a^2 – x^2) + (a^2 – x^2 )/2 ] dx= 2ª^3 / 3 Escríbela con símbolos y no es tan complicada

¡Ojalá lo hayas entendido!

Answer 7

Suena como si fuera el fin del mundo, pero en realidad es muy facil, que no te asuste que la integral es doble, es muy facil de resolver, solo es cuestion de resolver primero la integral de mas adentro y la respuesta que obtengas de esta, la integras por la variable que tenias de la integral de mas afuera ok
Un ejemplo a ver! Uno facil
SS2xdxdy
Primero integras con respecto a x, y lo que te salgas despues lo integras con respecto a y, asi de facil.Si? Solo imagina que las S son de esas integrales larguitas ok

Answer 8

busca aqui
http://huitoto.udea.edu.co/Matematicas/4.11.html