18 votos

Polinomios de ondulación

Me gustaría ser capaz de construir polinomios $p$ cuyas gráficas se parecen a este:

enter image description here

Podemos suponer que el intervalo de interés es $[-1, 1]$. Los requisitos en $p$ son:

(1) Equi-oscilación (o aproximadamente igual, de todos modos) entre los dos extremos. Una variación de 10% en los valores de los extremos estaría bien.

(2) los valores Cero y derivados en los extremos del intervalo, es decir, $p(-1) = p(1) =p'(-1) = p'(1) = 0$

Yo quiero hacer esto por grados hasta alrededor de los 30 o así. Sólo hasta grados sería ACEPTAR.

Si ayuda, estas cosas son un poco como los polinomios de Chebyshev (pero diferentes en los extremos).

El de la foto, tiene por ecuación $0.00086992073067855669451 - 0.056750328789339152999 t^2 + 0.60002383910750621904 t^4 - 2.3217878459074773378 t^6 + 4.0661558859963998471 t^8 - 3.288511471137768132 t^{10} + t^{12}$

Me llegó a través de la fuerza bruta métodos numéricos (resolver un sistema de ecuaciones no lineales, y después de hacer un montón de trabajo para encontrar buenos puntos de partida para la iteración). Estoy buscando un enfoque más inteligente y más fácil de implementar en el código.

Aquí es una idea que podría funcionar. Supongamos que queremos un polinomio de grado $n$. Comience con el polinomio de Chebyshev $T_{n-2}(x)$. Deje $Q(x) = T_{n-2}(sx)$, donde el factor de escala $s$ es elegido de manera que $Q(-1) = Q(1) = 0$. A continuación, vamos a $R(x) = (1-x^2)Q(x)$. Esto cumple con todos los requisitos, excepto que sus oscilaciones son muy desiguales, son muy pequeñas cerca de $\pm1$ demasiado grandes y cerca de cero. Redistribuir las raíces de $R$ un poco (de alguna manera??) a nivel de las oscilaciones.

Comentarios sobre respuestas

El uso de la técnica sugerida por achille hui en una respuesta de abajo, podemos muy fácilmente construir un polinomio con la forma deseada. Aquí hay uno:

achille hui solution

El único problema es que yo estaba esperando un polinomio de grado 12, y este tiene grado 30.

También, yo estaba esperando que la solución a crecer monótonamente fuera del intervalo de $[-1,1]$, y esto no se, como se puede ver aquí:

behaviour beyond unit interval

7voto

Joe Gauterin Puntos 9526

A partir de cualquier polinomio de Chebyshev de % bueno $1^{st}$% #% y cualquier raíz positiva $T_n(x), n > 2$. Composición $a$ con cualquier polinomio $T_n(ax)$ que es monótono en $q(x)$ $[-1,1]$ que y $q(-1) = -1, q(1) = 1$, entonces el $q'(-1) = q'(1) = 0$ es algo que usted quiere. Por ejemplo, $$ T_n\left(\cos\left(\frac{\pi}{2n}\right) \frac{x(3-x^2)}{2}\right) $$

3voto

alberta Puntos 16

Aquí el Newton de la versión de la Asíntota de código que realiza el trabajo. Ahora corre muy rápido en el rango solicitado. El polinomio es codificada por sus raíces (almacenado en la x de la matriz). Gracias al alma caritativa que tuvo la paciencia de tipo 4 espacios en el principio de cada línea :).

Edit. Desde bubba quería saber donde se encuentran los principales iteración paso vino, aquí una breve explicación. Primero de todo, como he dicho, es conveniente trabajar con la suma de $f(x)=\sum_k\log|x-x_k|$, en lugar del producto $\prod_k (x-x_k)$ (por tanto numérica y analítica razones). Esta suma inmersiones a $-\infty$ a cada raíz $x_k$ y alcanza un máximo local $y_{k-1}$ $y_k$ en algunos puntos de $z_{k-1}\in(x_{k-1},x_k)$$z_k\in(x_k,x_{k+1})$. Ahora, vamos a hacer un salto de fe y asumir que $y_{k-1}$ $y_k$ no son muy diferentes y también que el $z$-los puntos son aproximadamente la mitad de los intervalos correspondientes (la primera suposición de que va a suceder tarde o temprano, si nuestro algoritmo funciona en absoluto, y la segunda es en realidad bien justificada por sí misma, sino más razonable que cualquier otra conjetura acerca de la ubicación). Queremos cambiar la raíz de $x_k$, de modo que la maxima $y_{k-1}$ $y_k$ va a ser igual. Alinear y suponiendo que los intervalos de $(x_{k-1},x_k)$ $(x_k,x_{k+1})$ son de aproximadamente la misma longitud, vemos que las derivadas parciales de $f(z_{k-1})$ $f(z_k)$ con respecto al $x_k$ trata $\frac{1}{4(x_{k+1}-x_{k-1})}$ $-\frac{1}{4(x_{k+1}-x_{k-1})}$ respectivamente, por lo que, la solución de la correspondiente ecuación lineal, se obtiene la corrección de $p=(y_k-y_{k-1})(x_{k+1}-x_{k-1})/8$ a aplicar a $x_k$. Que explica la principal fórmula. Por desgracia, si simplemente usamos es, literalmente, un montón de cosas locas que puede suceder hasta la pérdida de la orden de las raíces, por lo que debemos asegurarnos de que nuestros cambios no son demasiado grandes en las primeras etapas, cuando los desequilibrios son bastante grandes. La segunda línea es una mera truncamiento de los turnos para el tamaño de lo que, creemos, se constituyen sólo una fracción de la distancia entre las raíces (que inicialmente es de $2/n$), por lo que no reordenamiento de las raíces o el atornillado de la iteración de Newton esquema se va a producir. La elección exacta de la constante de $0.3$ es empírica (es decir, acabo de comprobar que funcionaba en el rango solicitado y, digamos, $0.5$ no). La justificación formal del algoritmo requieren de una cuidadosa estimaciones (usted puede notar que algunos de los supuestos que he hecho están en la frontera de la ilusión), pero, desde todo lo que necesitamos es una secuencia de polinomios de no demasiado alto grado, pensé que valdría la pena intentar, sin preocuparse demasiado acerca de pruebas ya que el resultado final es verificable y una vez que lo tienes, a quién le importa qué medidas se han llevado a ello? Yo sé que esto es una triste cosa es decir para un matemático, pero, ya que estoy usando un programador de sombrero de aquí, Decidí que podría tratar de salirse con la suya :).

int n=51;
//The number of intermediate roots plus 1 (degree-3)

real[] x,y,z;
//The root array, the maximum array. and the critical point array

for(int k=0;k<n+1;++k) {x[k]=2*k/n-1; if(k>0) z[k-1]=(x[k-1]+x[k])/2;} 
//Just initialized the roots to an arithmetic progression
//and the critical points to midpoints between roots

real f(real t)
{
real s=2*log(1-t^2);
for(int k=1;k<n;++k) s+=log(abs(t-x[k]));
return s;
}
//This is just the logarithm of the absolute value of the polynomial
//with double roots at +1,-1 and simple roots at x[k], k=1,...,n-1

real g(real t)
{
real s=2*(1/(t+1)+1/(t-1));
for(int k=1;k<n;++k) 
s+=1/(t-x[k]);
return s;
}
//This is the derivative of f

real G(real t)
{
real s=-2*(1/(t+1)^2+1/(t-1)^2);
for(int k=1;k<n;++k) 
s-=1/(t-x[k])^2;
return s;
}
//This is the derivative of g

for(int m=0;m<15*n+30;++m)
//The outer cycle; the number of iterations is sufficient 
//to cover n up to 70 with the roots found with machine precision (1E-15) 
{
  for(int k=0;k<n/2;++k) 
  {
    real a=z[k];
    a-=g(a)/G(a); 
    y[k]=f(a); y[n-1-k]=y[k];
    z[k]=a; z[n-1-k]=-a;
  }
  //Newton update of critical points with symmetry taken into account

  real[] xx=copy(x);
  for (int k=1;k<n/2;++k) 
  { 
    real p=(y[k]-y[k-1])*(xx[k+1]-xx[k-1])/8; 
    if (abs(p)>0.3/n) p*=0.3/n/abs(p);
    x[k]+=p; x[n-k]-=p;   
  }
//The main iteration step: move each roots to balance the maxima
//adjacent to it. It can be done 
//better if we look beyond the adjacent maxima to evaluate what will 
//happen but, since it works in a reasonable time, I was too lazy to bother.

}

write(x);
write(y);
//outputs the roots and the maxima of f. Note that the roots at +1 and -1
//are double and the rest are simple.

pause();
//Just doesn't let the window close before you look at it.

i-Ciencias.com

I-Ciencias es una comunidad de estudiantes y amantes de la ciencia en la que puedes resolver tus problemas y dudas.
Puedes consultar las preguntas de otros usuarios, hacer tus propias preguntas o resolver las de los demás.

Powered by:

X