{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# MAP2220 - Fundamentos de Análise Numérica\n",
"## 2º Semestre de 2020 - Prof. Nelson Kuhl\n",
"## Exemplo com relação de recorrência em domínio discreto"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Iremos demonstrar o uso de relações de recorrência com polinômios ortogonais mônicos em domínios discretos. O produto interno é da forma $\\langle u, v\\rangle = \\sum_{i=1}^n \\omega_iu(x_i)v(x_i)$. Os polinômios até grau $m$,\n",
"onde $m \\le n-1$, são gerados por\n",
"$$\n",
" p_{k+1}(x) = (x - \\alpha_k)p_k(x) - \\beta_kp_{k-1}(x), \\quad 0 \\le k \\le m -1,\n",
"$$\n",
"com $p_{-1}(x) = 0$, $p_0(x) = 1$ e\n",
"$$\n",
" \\alpha_k = \\frac{\\langle xp_k,p_k\\rangle}{\\langle p_k, p_k\\rangle}, \\quad 0 \\le k \\le m-1, \\quad\n",
" \\beta_k = \\frac{\\langle p_k, p_k\\rangle}{\\langle p_{k-1}, p_{k-1}\\rangle}, \\quad 1 \\le k \\le m-1.\n",
"$$\n",
"Como são necessários apenas os valores dos polinômios nos pontos $x_i$, $1 \\le i \\le n$, somente estes são calculados."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"monicos (generic function with 1 method)"
]
},
"execution_count": 1,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"function monicos(m, xi, wi)\n",
" #= Dados os pontos x1, ..., xn e respectivos pesos w1, ..., wn,\n",
" o algoritmo \"contrói\" os m+1 (i.e, até grau m) primeiros\n",
" polinômios ortogonais em relação ao produto interno discreto\n",
" com peso. São retornados: os valores dos polinômios em uma\n",
" matriz n x (m+1), onde a coluna j refere-se ao grau (j-1);\n",
" as normas ao quadrado dos polinômios e os coeficientes da\n",
" relação de recorrência. =#\n",
"\n",
" n = length(xi); # quantidade de pontos\n",
" \n",
" if m > n - 1\n",
" print(\"O valor de m deve ser menor ou igual a $(n - 1)\")\n",
" return (nothing, nothing, nothing, nothing)\n",
" end\n",
"\n",
" PK = zeros(n, m+1) # para os polinômios\n",
" pkpk = zeros(m + 1); α = zeros(m); β = zeros(m) # para as normas e os coeficientes\n",
" \n",
" PK[:,1] = ones(n); pkpk[1] = sum(wi) # grau zero\n",
" \n",
" α[1] = sum(wi .* xi)/pkpk[1]; PK[:,2] = xi .- α[1]; pkpk[2] = sum(wi .* PK[:,2] .* PK[:,2]) # grau 1\n",
" \n",
" for k in 2:m # Graus 2 até m+1\n",
" α[k] = sum(wi .* xi .* PK[:,k] .* PK[:,k])/pkpk[k]\n",
" β[k] = pkpk[k]/pkpk[k-1]\n",
" PK[:,k+1] = (xi .- α[k]) .* PK[:,k] - β[k] * PK[:,k-1]\n",
" pkpk[k+1] = sum(wi .* PK[:,k+1] .* PK[:,k+1])\n",
" end\n",
" return (PK, pkpk, α, β)\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Uma vez gerados estes polinômios (na verdade, os valores destes polinômios nos pontos da malha), podemos usá-los para aproximar uma função $f$, tabelada nos pontos da malha, por um polinômio de grau menor ou igual a $m$, pelo MMQ. A solução é $g(x) = \\sum_{k=0}^m c_kp_k(x)$ com\n",
"$$\n",
" c_k = \\frac{\\langle p_k, f\\rangle}{\\langle p_k, p_k\\rangle} =\n",
" \\frac{\\sum_{i=1}^n \\omega_ip_k(x_i)f(x_i)}{\\sum_{i=1}^n \\omega_ip_k^2(x_i)}, \\quad 0 \\le k \\le m.\n",
"$$\n",
"Note que para os cálculos dos coeficientes, **usamos somente** os valores dos polinômios nos pontos da tabela e as suas normas, que são calculados na função *monicos* acima. Se quisermos calcular o valor de $g(x)$ em um ponto $x$, pdemos fazê-lo sen conhecer as fórmulas explícitas dos $p_k(x)$, usando o algoritmo de Clenshaw. Tudo o que precisamos são os coeficientes $c_k$ e os valores de $\\alpha_k$ e $\\beta_k$, já calculados."
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"clenshaw (generic function with 1 method)"
]
},
"execution_count": 2,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Algoritmo de Clenshaw adaptado para a recursão dos polinômios mônicos\n",
"function clenshaw(x, c, alfa, beta)\n",
" a = x .- alfa\n",
" M = length(c)\n",
" y = zeros(M)\n",
" y[M] = c[M]; y[M-1] = c[M-1] + a[M-1]*y[M]\n",
" for k in M-2:-1:1\n",
" y[k] = c[k] + a[k]*y[k+1] - beta[k+1]*y[k+2]\n",
" end\n",
" return y[1]\n",
"end"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Exemplo\n",
"Para praticarmos os métodos acima, faremos alguns exemplos arbitrários."
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [],
"source": [
"using LinearAlgebra, Plots"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×9 Array{Float64,2}:\n",
" 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9\n",
" 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0"
]
},
"execution_count": 4,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Pontos e pesos\n",
"x = collect(0.1:0.1:0.9)\n",
"w = ones(length(x)) # Pesos iguais a 1\n",
"[x'; w']"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [],
"source": [
"# Construção da matriz dos polinômios, respectivas normas ao quadrado e coeficientes da recorrência\n",
"(MP, normp2, α, β) = monicos(4, x, w); # Grau menor ou igual a 4"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9×5 Array{Float64,2}:\n",
" 1.0 -0.4 0.0933333 -0.0168 0.0024\n",
" 1.0 -0.3 0.0233333 0.0084 -0.0036\n",
" 1.0 -0.2 -0.0266667 0.0156 -0.00188571\n",
" 1.0 -0.1 -0.0566667 0.0108 0.00154286\n",
" 1.0 0.0 -0.0666667 7.40149e-18 0.00308571\n",
" 1.0 0.1 -0.0566667 -0.0108 0.00154286\n",
" 1.0 0.2 -0.0266667 -0.0156 -0.00188571\n",
" 1.0 0.3 0.0233333 -0.0084 -0.0036\n",
" 1.0 0.4 0.0933333 0.0168 0.0024"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"MP # valores dos polinômios nos pontos; a coluna k é referente a p_{k-1}"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 7,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Verificação aproximada da ortogonalidade\n",
"isapprox(MP'*(diagm(w)*MP), diagm(normp2))"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×4 Array{Float64,2}:\n",
" 0.5 0.5 0.5 0.5\n",
" 0.0 0.0666667 0.0513333 0.0462857"
]
},
"execution_count": 8,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[α'; β'] # Parâmetros da recorrência"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 9,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# f gerada aleatoriamente segundo a distribuição normal\n",
"f = randn(length(x))\n",
"scatter(x, f, label = \"dados\", legend = :outertopright)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×5 Adjoint{Float64,Array{Float64,1}}:\n",
" 0.0703584 -0.316237 -1.09978 16.9108 -87.3616"
]
},
"execution_count": 10,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Cálculo dos coeficientes c_k\n",
"c = zeros(length(normp2))\n",
"for k in 1:length(normp2)\n",
" c[k] = sum(w .* MP[:,k] .* f)/normp2[k]\n",
"end\n",
"c'"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {},
"outputs": [],
"source": [
"# Calculo de g em vários pontos para obter o seu gráfico\n",
"xx = collect(range(x[1], stop=x[end], step=0.01))\n",
"N = length(xx)\n",
"g = zeros(N)\n",
"for i in 1:N\n",
" g[i] = clenshaw(xx[i], c, α, β)\n",
"end"
]
},
{
"cell_type": "code",
"execution_count": 12,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 12,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"plot!(xx, g, label = \"aprox\")"
]
},
{
"cell_type": "code",
"execution_count": 13,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9-element Array{Float64,1}:\n",
" 1.0\n",
" 0.5\n",
" 0.3333333333333333\n",
" 0.25\n",
" 0.2\n",
" 0.16666666666666666\n",
" 0.14285714285714285\n",
" 0.125\n",
" 0.1111111111111111"
]
},
"execution_count": 13,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Agora usando pesos, apenas por curiosidade.\n",
"w = 1 ./ collect(1:length(x)) # w_i = 1/i"
]
},
{
"cell_type": "code",
"execution_count": 14,
"metadata": {},
"outputs": [],
"source": [
"# Construção da matriz dos polinômios, respectivas normas ao quadrado e coeficientes da recorrência\n",
"(MP, normp2, α, β) = monicos(4, x, w); # Grau menor ou igual a 4"
]
},
{
"cell_type": "code",
"execution_count": 15,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"9×5 Array{Float64,2}:\n",
" 1.0 -0.218137 0.039964 -0.00604243 0.000746506\n",
" 1.0 -0.118137 -0.0166937 0.0110894 -0.00277325\n",
" 1.0 -0.0181372 -0.0533513 0.0125264 -0.000350327\n",
" 1.0 0.0818628 -0.070009 0.00426862 0.00260582\n",
" 1.0 0.181863 -0.0666667 -0.00768398 0.00308571\n",
" 1.0 0.281863 -0.0433243 -0.0173314 0.000479897\n",
" 1.0 0.381863 1.79972e-5 -0.0186736 -0.0034211\n",
" 1.0 0.481863 0.0633603 -0.00571061 -0.00442675\n",
" 1.0 0.581863 0.146703 0.0275576 0.00405349"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"MP # valores dos polinômios nos pontos; a coluna k é referente a p_{k-1}; observe a diferença"
]
},
{
"cell_type": "code",
"execution_count": 16,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"true"
]
},
"execution_count": 16,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Verificação aproximada da ortogonalidade\n",
"isapprox(MP'*(diagm(w)*MP), diagm(normp2))"
]
},
{
"cell_type": "code",
"execution_count": 17,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"2×4 Array{Float64,2}:\n",
" 0.318137 0.548439 0.518164 0.516837\n",
" 0.0 0.0578573 0.0489099 0.044345"
]
},
"execution_count": 17,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"[α'; β'] # Parâmetros da recorrência; observe a diferença"
]
},
{
"cell_type": "code",
"execution_count": 18,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"1×5 Adjoint{Float64,Array{Float64,1}}:\n",
" 0.0235932 0.257146 -4.29747 27.7433 -110.062"
]
},
"execution_count": 18,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Cálculo dos coeficientes c_k; observe a diferença\n",
"c = zeros(length(normp2))\n",
"for k in 1:length(normp2)\n",
" c[k] = sum(w .* MP[:,k] .* f)/normp2[k]\n",
"end\n",
"c'"
]
},
{
"cell_type": "code",
"execution_count": 19,
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n"
]
},
"execution_count": 19,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# Cálculo da nova aproximação nos pontos xx anteriores e respectivo gráfico\n",
"gw = zeros(N)\n",
"for i in 1:N\n",
" gw[i] = clenshaw(xx[i], c, α, β)\n",
"end\n",
"plot!(xx, gw, label = \"aprox c/ w\")"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Julia 1.5.2",
"language": "julia",
"name": "julia-1.5"
},
"language_info": {
"file_extension": ".jl",
"mimetype": "application/julia",
"name": "julia",
"version": "1.5.2"
}
},
"nbformat": 4,
"nbformat_minor": 2
}