Problemas de contorno en una dimensión

Estudiamos problemas de contorno para ecuaciones diferenciales ordinarias. La variable respecto de la que derivamos típicamente no representa el tiempo, sino que modelizamos problemas estacionarios en una dimensión espacial. Aunque se pueden definir problemas de contorno de interés con la variable temporal

Comenzamos con una sóla ecuación diferencial: $$ y'' = f(x,y,y') $$ Definida en un intervalo $[a,b]$, precribiendo los valores de $f$ en los extremos del intervalo: $$ y'' = f(x,y,y'),\; a\leq x\leq b, y(a)=\alpha, y(b)=\beta $$

Comparamos nuestro problema de contorno (PC): $$ y'' = f(x,y,y'),\; a\leq x\leq b, y(a)=\alpha, y(b)=\beta $$ con el problema de valor de inicial (PVI) para la misma ecuación diferencial, que necesita que especifiquemos el valor inicial $y(a)$ y la velocidad inicial $y'(a)$: $$ y'' = f(x,y,y'),\; a\leq x\leq b, y(a)=\alpha, y'(a)=t $$

El problema de valor inicial tiene solución única si la derivada $f_y$ de $f$ respecto de su segundo argumento ($y$), y la derivada $f_{y'}$ respecto de su tercer argumento ($y'$), están acotadas.

El problema de contorno necesita condiciones adicionales, por ejemplo que la derivada $f_{y'}$ sea positiva.

Método del disparo

Comenzamos con esta observación:

Si la solución de (PC) es $y$, y elegimos $t=y'(a)$ entonces $y$ también es la única solución de (PVI).

Por tanto, si encontramos $t$, hemos resuelto nuestro problema de contorno (PC).

Definimos una función $D(t)$:

$$ D(t) = y(b),\; y \text{ es la única solución de (PVI) con } y'(a)=t $$

La solución $y$ a (PC) es la única solución al (PVI) con $y(a)=\alpha, y'(a)=t$ tal que $D(t)=\beta$.

La función $D$ es una función que puede ser no lineal, y tenemos que encontrar una raíz de $D$, con cualquiera de los métodos que estudiamos al principio de la asignatura.

Ejemplo

Representamos el desplazamiento vertical de una viga horizontal mediante una función $w(x)$ que indica el desplazamiento vertical en la posición $x$ medida desde el extremo izquierdo de la viga.

viga

La siguiente ecuación modeliza el desplazamiento vertical $w$ de una viga con carga homogénea y sujeta en los dos extremos:

$$ w''(x) = \frac{S}{EI}w(x) + \frac{q}{2EI}x(x-l)\; 0<x<l $$

donde

  • $l$: longitud de la viga
  • $E$: módulo elástico de la viga
  • $I$: momento de inercia de la sección de la viga
  • $S$: tensión horizontal en los extremos de la viga
  • $q$: carga de la viga por unidad de longitud

Fijamos los valores $l=120, E=3\cdot 10^{7}, S=1000, I=625$ (en unidades compatibles que no especificamos).

Comenzamos por definir la función fun que necesitamos para resolver los (PVI)

In [4]:
l = 120;
Q = 100;
E = 3e7;
S = 1e3;
I = 625;

fun = @(x,ys)([ys(2), ((S/(E*I))*ys(1)  + (Q*x/(2*E*I))*(x-l))]);

Definimos una función disparo, que recibe como argumento la velocidad del (PVI), y devuelve $D(t)-\beta$, de modo que la velocidad t que buscamos para el (PVI) es una raíz de la función disparo que podemos buscar con root, o cualquier otro método (bisección, newton, secante, etc).

In [5]:
%Esta linea define x0, xf (numeros) y x_span (una lista con esos dos mismos numeros)
x0 = 0;
xf = l;
x_span = [x0, xf];
alpha = 0;
beta = 0;
In [6]:
function target = NLS(v)
    l = 120;
    Q = 100;
    E = 3e7;
    S = 1e3;
    I = 625;
    x0 = 0;
    xf = l;
    funNL = @(x,ys)([
        ys(2), 
        ((S/(E*I))*ys(1) + (Q/(2*E*I))*x*(x-l))*(1 + ys(2)^2)^1.5
    ]);
    x_span = [0, l];
    alfa = 0;
    beta = 0;
    y0 = [alfa, v];
    [t,y] = ode45(funNL, x_span, y0);
    target = y(length(t), 1) - beta;
end
In [7]:
vel_x0 = fsolve(@NLS, 0)

vel_x0 =    3.8397e-04

Ahora recogemos la solución completa del (PVI) con esa velocidad t que acabamos de encontrar.

In [8]:
l = 120;
Q = 100;
E = 3e7;
S = 1e3;
I = 625;
funNL = @(x,ys)([
    ys(2), 
    ((S/(E*I))*ys(1) + (Q/(2*E*I))*x*(x-l))*(1 + ys(2)^2)^1.5
    ]);

nx = 200;
xs = linspace(x0,xf,nx);
y0 = [alpha, vel_x0];
[t,ys] = ode45(funNL, xs, y0);
ys_NL = ys(:,1);
In [9]:
hold on;
plot(xs, -transpose(ys_NL), 'g');
xlabel('x');
ylabel('desplazamiento vertical');
title('Solucion del problema de contorno');
hold off;

Ejercicio

La ecuación de van der Pol modeliza ciertos circuitos eléctricos, entre otras cosas:

The Van der Pol oscillator was originally proposed by the Dutch electrical engineer and physicist Balthasar van der Pol while he was working at Philips. Van der Pol found stable oscillations, which he subsequently called relaxation-oscillations and are now known as a type of limit cycle in electrical circuits employing vacuum tubes.

La ecuación para la corriente que fluye por cierta parte del circuito es: $$ {\displaystyle {d^{2}x \over dt^{2}}-\mu (1-x^{2}){dx \over dt}+x=0.} $$

  • Para $\mu = 8.5$, $A = 1.2$, $\omega = \pi/5$, resolved el problema de valor inicial $x(0)=0, x'(0)=1$, hasta alcanzar $t=50$.
  • Para $\mu=0.5$, resolved el problema de contorno $x(0)=0,x(2)=1$.
In [ ]:
In [ ]:

Ejercicio

Añadimos una fuente de corriente alterna al circuito, de modo que ahora la ecuación es: $$ {\displaystyle {d^{2}x \over dt^{2}}-\mu (1-x^{2}){dx \over dt}+x=A\sin(\omega t).} $$

  • Para $\mu = 8.5$, $A = 1.2$, $\omega = \pi/5$, resolved el problema de valor inicial $x(0)=0, x'(0)=1$, hasta alcanzar $t=50$.
  • Para $\mu=0.5$, resolved el problema de contorno $x(0)=0,x(2)=1$.
In [ ]:
In [ ]:

Ejercicio

Resolved el mismo problema de contorno del problema anterior, pero usando el método bvp4c de matlab.

In [3]:
In [ ]:
In [ ]: