Utilización del método de colocación para resolver una EDO no lineal de valores de contorno

Question

Tengo la siguiente EDO

$$ u'' = -(1 + e^u), \quad u(0)=0,\quad u(1)=1$$

Quiero dividir el intervalo $[0,1]$ en $n-1$ subintervalos iguales, cada uno con una longitud $h=1/(n-1)$ y tomamos la solución aproximada

$$ v(t,x) = \sum_{i=1}^{n} x_i t^{i-1} $$

polinomio de grado $n-1$ así que $v(0,x)=0$ y $v(1,x)=1$ (deben cumplirse las condiciones de contorno).

Por lo tanto, nuestro objetivo es determinar $x_2,x_3,...,x_{n-1}$ puntos interiores puesto que ya sabemos $x_1 = 0$ . Una vez encontrados, tenemos una aproximación a nuestra solución verdadera $u(t) \approx v(t,x)$ .

Mi pregunta es, ¿cómo podemos implementar esto en matlab? porque resolver tal sistema es bastante tiempo a mano.

Answer 1

Para realizar Colocación introducimos la discretización regular $t_i = i/n$ con $0\leq i\leq n$ del dominio $[0,1]$ (es decir, un conjunto de $n+1$ puntos de colocación). A continuación, aproximamos $u$ por la función polinómica $p(t) = \sum a_k t^k$ de grado $n$ que satisface el problema anterior en los puntos de colocación. Por lo tanto, $-p''(t_i) = 1 + \exp\!\big(p(t_i)\big)$ en el dominio interior $1\leq i \leq n-1$ , y $p(t_0)=0$ , $p(t_n) = 1$ en los límites. En cuanto a los coeficientes $a_k$ este sistema se escribe como $$ -\sum_{k=0}^{n-2} (k+2) (k+1) a_{k+2} \frac{i^k}{n^k} = 1 + \exp\left(\sum_{k=0}^n a_k \frac{i^k}{n^k}\right) , \qquad 1\leq i \leq n-1 $$ y $a_0 = 0$ , $\sum a_k = 1$ . Este sistema algebraico no lineal puede entonces resolverse numéricamente en términos de los coeficientes $a_k$ por ejemplo, utilizando el programa de Matlab fsolve .

Por ejemplo, si $n=2$ consideramos los puntos de colocación $\lbrace 0, 0.5, 1\rbrace$ y buscamos los coeficientes $a_0$ , $a_1$ , $a_2$ . Estos coeficientes satisfacen $a_0 = 0$ y el sistema algebraico $$ \left\lbrace \begin{aligned} &{-2} a_{2} = 1 + \exp\left(\tfrac{1}{2} a_1 + \tfrac{1}{4} a_2\right)\\ &a_1 + a_2 = 1 \end{aligned}\right. $$ cuyas soluciones reales se obtienen numéricamente $$ (a_1,a_2) \in \lbrace (2.78945, -1.78945), (10.6101, -9.61014) \rbrace . $$ se representan a continuación:

Se observa que estas soluciones no son necesariamente únicas. Comparación con Matlab bvp4c se da a continuación, donde el algoritmo se inicializa con la constante $0.5$ :

sol = bvp4c(@(x,y) [y(2) -(1+exp(y(1)))], @(ya,yb) [ya(1) yb(1)-1], bvpinit(linspace(0,1,10), [0.5; 0.5]));
plot(sol.x,sol.y(1,:));
xlabel('t');
ylabel('u');

Si la inicialización se realiza con la constante $2.5$ , entonces se obtiene la solución superior.

Answer 2

El polinomio se puede determinar, por ejemplo, con el código (en octava, puede haber ligeras diferencias para matlab).

• En primer lugar, define una función que toma una matriz de coeficientes y devuelve una matriz que contiene las condiciones de contorno y los residuos de la EDO en los puntos de colocación.

  function f = system(c)
    % degree of the interpolating polynomial
    % is also the number of subdivisions
    deg = length(c)-1;
    f = zeros(1,deg+1);
    % left boundary condition
    f(1)=polyval(c,0);
    % second derivative of polynomial
    c2 = polyder(polyder(c));
    % enforce differential equation at inner points
    for k = 1:deg-1
      xk = k*1.0/deg;
      f(k+1)= polyval(c2,xk) + (1+exp(polyval(c,xk)));
    end
    % right boundary condition
    f(deg+1)=polyval(c,1) - 1;
  end

• Esta función del sistema no lineal se puede llamar ahora en fsolve donde la matriz de inicialización también determina el grado del polinomio de colocación

  clf;
  % define target accuracy
  options = optimset("TolX", 1e-9, "TolFun", 1e-6);
  % number of subdivisions = polynomial degree
  deg = 5;
  hold on;
  % loop over different initializations to find multiple solutions
  for k = 0:5
     c = zeros(1,deg+1);
     c(deg) = 1+3*k; c(deg-1) = -3*k;
     % log the initial coefficients
     c
     % solve system and log the result
     c = fsolve(@(c) system(c), c, options)
     % plot the solution
     x = linspace(0,1,101);
     plot(x, polyval(c,x))
  end
  hold off;

La ejecución de este código encuentra dos soluciones, que se repiten 3 veces para las inicializaciones elegidas:

 0.24661  -0.38053    -0.31725   -1.02286    2.47402    0.00000
30.41419  -49.23356   11.70577   -1.37713    9.49074    0.00000

con la trama

---

Adenda: Tomando las soluciones de grado 2 como en la respuesta de Harry49

 -1.78945    2.78945    0.00000
 -9.61015   10.61015    0.00000

y utilizando estos polinomios y su derivada para inicializar el solucionador BVP (en python, ya que mi versión de octave no tiene esta funcionalidad), da el siguiente gráfico, donde los polinomios se dibujan en gris debajo de las soluciones más exactas:

Answer 3

Existe una solución analítica de la EDO en cuadraturas.

Dejemos que $P(u)=\dfrac{du}{dx}.$

Entonces $$u''=\dfrac{dP}{dt} = \dfrac{dP}{du}\dfrac{du}{dx} = PP',$$ $$P^2 = -2(u+e^u+\mathrm{const}),$$ $$u'^2 = 2(C_1-u-e^u),$$ $$u' = \pm\sqrt2\sqrt{C_1-u-e^u},\tag1$$ $$x=C_2\pm\dfrac1{\sqrt2}\int \dfrac{du}{\sqrt{C_1 -u - e^u}}.\tag2$$ Asumiendo la derivada $P$ positivo y teniendo en cuenta las condiciones de contorno, la ecuación $(2)$ puede presentarse en forma de $$x=I(u,C_1)=\dfrac1{\sqrt2}\int\limits_0^u \dfrac{dv}{\sqrt{C_1 -v - e^v}},\tag3$$ donde la constante $C_1$ puede definirse a partir de la ecuación $$I(1,C_1) = 1.\tag4$$ Teniendo en cuenta que $C_1\ge e+1$ y utilizando la sustitución $w=1-v,$ se puede conseguir $$I(1,e) = \dfrac1{\sqrt2}\int\limits_0^1\dfrac{dv}{\sqrt{1-v+e(1-e^{v-1})}} =\dfrac1{\sqrt2}\int\limits_0^1\dfrac{dw}{\sqrt{w +e(1-e^{-w})}}$$ $$\le \dfrac1{\sqrt2}\int\limits_0^1\dfrac{dw}{\sqrt{w +e(w-\frac12w^2)}} = 2\int\limits_0^1\dfrac{d\sqrt w}{\sqrt{2+2e-e w}} = \dfrac2{\sqrt e}\int\limits_0^1\dfrac{dz}{\sqrt{2+2e^{-1}-z^2}}$$ $$= \dfrac2{\sqrt e}\arcsin\dfrac1{\sqrt{2+2e^{-1}}}]\approx 0.788<1$$ ( Cálculos numéricos de Wolfram Alpha)dv) dar $I(1,e)\approx 0.776236$ )

Este resultado se contradice con $(4),$ por lo que las condiciones de contorno $u(0)=0$ y $u(1)=1\quad \color{green}{\mathbf{cannot\,be\,satisfied\,in\,the\,model\, with\,the\,constant\,sign\,of\, P}}$ (ver también los comentarios de $\color{green}{\textbf{LutzL}}).$

Esto significa que la función $u(x)$ tiene el máximo global, que coincide con el punto de ramificación de $x(u).$

Dejemos que $(x_m,u_m)$ es el punto de ramificación, en el que $$u'=\begin{cases} \sqrt2\sqrt{C_1-u-e^u},\quad\text{if}\quad u < u_m\\[4pt] 0,\quad\text{if}\quad u = u_m\\[4pt] -\sqrt2\sqrt{C_1-u-e^u},\quad\text{if}\quad u > u_m. \end{cases}$$

$$u_m>1,\quad\dfrac{dx}{du}\bigg|_{\large u_m} =\infty,\quad C_1 = x_m = u_m+e^{u_m},\quad u_m = x_m-W(e^{x_m}),\tag6$$

donde $W(x)$ es Función Lambert W .

Entonces $$x_{1,2}(u) = x_m \pm \int\limits_u^{u_m}\dfrac{dv}{\sqrt{2(x_m - v-e^v)}},\tag7$$

$$\int\limits_0^{u_m} \dfrac{dv}{\sqrt{2(x_m - v-e^v)}} + \int\limits_1^{u_m} \dfrac{dv}{\sqrt{2(x_m - v-e^v)}} = 1,$$ $$\int\limits_0^{1} \dfrac{dv}{\sqrt{2(x_m - v-e^v)}} + 2\int\limits_1^{u_m} \dfrac{dv}{\sqrt{2(x_m - v-e^v)}} = 1,\tag8$$ con el soluciones numéricas $$\dbinom{u_m}{x_m}\in\left\{\approx\dbinom{1.08770}{4.05514},\dbinom{3.67011}{42.92633}\right\}.\tag9$$