Baixe o app para aproveitar ainda mais
Prévia do material em texto
Computação Numérica Aula de Laboratório 5: Integração numérica Professor Filipe Taveiros Considere a seguinte função: Aquecimento Considere a seguinte função: 𝑓 𝑥 = 4𝑥3 ⋅ 𝑒−2𝑥 Aquecimento Considere a seguinte função: 𝑓 𝑥 = 4𝑥3 ⋅ 𝑒−2𝑥 É fácil de calcular a integral 𝑓(𝑥) 𝑑𝑥 ? Aquecimento Considere a seguinte função: 𝑓 𝑥 = 4𝑥3 ⋅ 𝑒−2𝑥 É fácil de calcular a integral 𝑓(𝑥) 𝑑𝑥 ? 𝑓(𝑥) 𝑑𝑥 = − 1 2 𝑒−2𝑥(3 + 6 𝑥 + 6 𝑥2 + 4 𝑥3) Aquecimento Considere a seguinte função: 𝑓 𝑥 = 4𝑥3 ⋅ 𝑒−2𝑥 É fácil de calcular a integral 𝑓(𝑥) 𝑑𝑥 ? 𝑓(𝑥) 𝑑𝑥 = − 1 2 𝑒−2𝑥(3 + 6 𝑥 + 6 𝑥2 + 4 𝑥3) Mais ou menos... E se não fosse? Aquecimento Aquecimento Vamos aquecer!! Aquecimento Vamos aquecer!! 1. Defina a função f(x) no Scilab Aquecimento Vamos aquecer!! 1. Defina a função f(x) no Scilab 2. Também defina a integral da função f(x) no Scilab. Vamos usá-la como referência. Aquecimento Vamos aquecer!! 1. Defina a função f(x) no Scilab 2. Também defina a integral da função f(x) no Scilab. Vamos usá-la como referência. 3. Calcule a integral de f(x) no intervalo [0,3], utilizando: a) A função analítica b) Uma função nativa do Scilab Aquecimento Aquecimento Scilab 5.5.1 Console Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); Scilab 5.5.1 Console Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); Scilab 5.5.1 Console Na hora de definir a função, nunca esquecer de utilizar o ponto “.” para as operações de multiplicação, divisão e exponenciação. Isto permite que a função receba um vetor como argumento de entrada!! Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); -->deff('y=intf(x)','y=-(1/2)*exp(-2*x).*(3+6*x+6*x.^2+4*x.^3)') Scilab 5.5.1 Console Na hora de definir a função, nunca esquecer de utilizar o ponto “.” para as operações de multiplicação, divisão e exponenciação. Isto permite que a função receba um vetor como argumento de entrada!! Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); -->deff('y=intf(x)','y=-(1/2)*exp(-2*x).*(3+6*x+6*x.^2+4*x.^3)') Scilab 5.5.1 Console Na hora de definir a função, nunca esquecer de utilizar o ponto “.” para as operações de multiplicação, divisão e exponenciação. Isto permite que a função receba um vetor como argumento de entrada!! Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); -->deff('y=intf(x)','y=-(1/2)*exp(-2*x).*(3+6*x+6*x.^2+4*x.^3)') Scilab 5.5.1 Console Na hora de definir a função, nunca esquecer de utilizar o ponto “.” para as operações de multiplicação, divisão e exponenciação. Isto permite que a função receba um vetor como argumento de entrada!! Função nativa do Scilab Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); -->deff('y=intf(x)','y=-(1/2)*exp(-2*x).*(3+6*x+6*x.^2+4*x.^3)') -->integrate('f(x)','x',0,3) ans = 1.2731942 Scilab 5.5.1 Console Na hora de definir a função, nunca esquecer de utilizar o ponto “.” para as operações de multiplicação, divisão e exponenciação. Isto permite que a função receba um vetor como argumento de entrada!! Função nativa do Scilab Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); -->deff('y=intf(x)','y=-(1/2)*exp(-2*x).*(3+6*x+6*x.^2+4*x.^3)') -->integrate('f(x)','x',0,3) ans = 1.2731942 Scilab 5.5.1 Console Na hora de definir a função, nunca esquecer de utilizar o ponto “.” para as operações de multiplicação, divisão e exponenciação. Isto permite que a função receba um vetor como argumento de entrada!! Função nativa do Scilab Teorema fundamental do cálculo Aquecimento -->deff('y=f(x)','y=(4*x.^3).*exp(-2*x)'); -->deff('y=intf(x)','y=-(1/2)*exp(-2*x).*(3+6*x+6*x.^2+4*x.^3)') -->integrate('f(x)','x',0,3) ans = 1.2731942 -->intf(3)-intf(0) ans = 1.2731942 Scilab 5.5.1 Console Na hora de definir a função, nunca esquecer de utilizar o ponto “.” para as operações de multiplicação, divisão e exponenciação. Isto permite que a função receba um vetor como argumento de entrada!! Função nativa do Scilab Teorema fundamental do cálculo Métodos estudados Métodos estudados 1. Trapézio Métodos estudados 1. Trapézio 2. 1/3 de Simpson Métodos estudados 1. Trapézio 2. 1/3 de Simpson 3. 3/8 de Simpson Trapézio Trapézio Scinotes Trapézio function i=trap(a,b) Scinotes Trapézio function i=trap(a,b) i=(b-a)*(f(b)+f(a))/2; Scinotes Trapézio function i=trap(a,b) i=(b-a)*(f(b)+f(a))/2; endfunction Scinotes Trapézio function i=trap(a,b) i=(b-a)*(f(b)+f(a))/2; endfunction Scinotes Todos devem copiar este código no Scilab! Trapézio Scilab 5.5.1 Console Trapézio -->trap(0,3) Scilab 5.5.1 Console Trapézio -->trap(0,3) ans = 0.4015579 Scilab 5.5.1 Console Trapézio -->trap(0,3) ans = 0.4015579 Scilab 5.5.1 Console O valor exato é 1.273! O que raio é isso?? Trapézio -->trap(0,3) ans = 0.4015579 Scilab 5.5.1 Console O valor exato é 1.273! O que raio é isso?? Trapézio Como melhorar? Trapézio Como melhorar? Segmentando o intervalo [𝑎, 𝑏] em N subintervalos! Trapézio Como melhorar? Segmentando o intervalo [𝑎, 𝑏] em N subintervalos! Para cada um dos subintervalos, calculamos a sua contribuição de área utilizando o método do trapézio. Trapézio Como melhorar? Trapézio 1. Façamos então um contador 𝑘 que vai de 0 até N-1. O espaço entre cada subintervalo é dado por ℎ = 𝑏 − 𝑎 𝑁 Trapézio 1. Façamos então um contador 𝑘 que vai de 0 até N-1. O espaço entre cada subintervalo é dado por ℎ = 𝑏 − 𝑎 𝑁 2. Assim, os limites dos subintervalos serão dados por 𝑎𝑘 = 𝑎 + ℎ𝑘 𝑏𝑘 = 𝑎 + ℎ(𝑘 + 1) Trapézio 1. Façamos então um contador 𝑘 que vai de 0 até N-1. O espaço entre cada subintervalo é dado por ℎ = 𝑏 − 𝑎 𝑁 2. Assim, os limites dos subintervalos serão dados por 𝑎𝑘 = 𝑎 + ℎ𝑘 𝑏𝑘 = 𝑎 + ℎ(𝑘 + 1) 3. Para cada subintervalo, calcule a sua área pela fórmula do trapézio. Trapézio 1. Façamos então um contador 𝑘 que vai de 0 até N-1. O espaço entre cada subintervalo é dado por ℎ = 𝑏 − 𝑎 𝑁 2. Assim, os limites dos subintervalos serão dados por 𝑎𝑘 = 𝑎 + ℎ𝑘 𝑏𝑘 = 𝑎 + ℎ(𝑘 + 1) 3. Para cada subintervalo, calcule a sua área pela fórmula do trapézio. 4. Some a área de todos os trapézios. Trapézio 1. Façamos então um contador 𝑘 que vai de 0 até N-1. O espaço entre cada subintervalo é dado por ℎ = 𝑏 − 𝑎 𝑁 2. Assim, os limites dos subintervalos serão dados por 𝑎𝑘 = 𝑎 + ℎ𝑘 𝑏𝑘 = 𝑎 + ℎ(𝑘 + 1) 3. Para cada subintervalo, calcule a sua área pela fórmula do trapézio. 4. Some a área de todos os trapézios. 5 minutos pra vocês fazerem sozinhos! Go, go, go!!! Trapézio Composto Trapézio Composto Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; i = 0; Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; i = 0; for k = 0:n-1 Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; i= 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (b-a)*(f(b)+f(a))/2; Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (b-a)*(f(b)+f(a))/2; end Scinotes Trapézio Composto function i=trapcomposto(a0, bf, n) h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (b-a)*(f(b)+f(a))/2; end endfunction Scinotes Trapézio composto Scilab 5.5.1 Console Trapézio composto Scilab 5.5.1 Console -->trapcomposto(0,3,2) Trapézio composto Scilab 5.5.1 Console -->trapcomposto(0,3,2) ans = 1.2089671 Trapézio composto Scilab 5.5.1 Console -->trapcomposto(0,3,2) ans = 1.2089671 Trapézio composto -->trapcomposto(0,3,3) ans = 1.2612942 Scilab 5.5.1 Console Trapézio composto -->trapcomposto(0,3,4) ans = 1.2664979 Scilab 5.5.1 Console Trapézio composto -->trapcomposto(0,3,5) ans = 1.2681641 Scilab 5.5.1 Console Trapézio composto -->trapcomposto(0,3,6) ans = 1.2692405 Scilab 5.5.1 Console Trapézio composto -->trapcomposto(0,3,7) ans = 1.270034 Scilab 5.5.1 Console Como melhorar? Como melhorar? IDEIA: Ao invés de interpolar dois pontos por uma reta, vamos interpolar por um polinômio de grau 2. Como melhorar? RELEMBRANDO: Para gerar um polinômio de grau 2, precisamos de 3 pontos. Portanto, usaremos os dois pontos do subintervalo mais um ponto extra no meio do subintervalo. Esta é a regra de 1/3 de Simpson. Como melhorar? RELEMBRANDO: Como melhorar? RELEMBRANDO: A fórmula do método 1/3 de Simpson é Como melhorar? RELEMBRANDO: A fórmula do método 1/3 de Simpson é 3 minutos pra vocês modificarem o método do trapézio composto para tornar-se o método 1/3 de Simpson. Go, go, go!!! 1/3 Simpson 1/3 Simpson Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; i = 0; Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; i = 0; for k = 0:n-1 Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (1/3)*((b-a)/2)*(f(a)+4*f((a+b)/2)+f(b)); Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (1/3)*((b-a)/2)*(f(a)+4*f((a+b)/2)+f(b)); end Scinotes 1/3 Simpson function i=simpson13(a0, bf, n) n=n/2; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (1/3)*((b-a)/2)*(f(a)+4*f((a+b)/2)+f(b)); end endfunction Scinotes 1/3 Simpson Scilab 5.5.1 Console 1/3 Simpson -->simpson13(0,3,2) Scilab 5.5.1 Console 1/3 Simpson -->simpson13(0,3,2) ans = 1.4781035 Scilab 5.5.1 Console 1/3 Simpson -->simpson13(0,3,2) ans = 1.4781035 Scilab 5.5.1 Console 1/3 Simpson -->simpson13(0,3,2) ans = 1.4781035 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos. 1/3 Simpson -->simpson13(0,3,4) ans = 1.2856748 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos. 1/3 Simpson -->simpson13(0,3,6) ans = 1.2718893 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos. 1/3 Simpson -->simpson13(0,3,8) ans = 1.2720081 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos. 1/3 Simpson -->simpson13(0,3,10) ans = 1.2725231 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos. Como melhorar? Como melhorar? IDEIA: Ao invés de interpolar dois pontos por um polinômio de grau 2, vamos interpolar utilizando um polinômio de grau 3. Como melhorar? RELEMBRANDO: Para gerar um polinômio de grau 3, precisamos de 4 pontos. Portanto, usaremos os dois pontos do subintervalo mais dois pontos extras. Esta é a regra de 3/8 de Simpson. Como melhorar? RELEMBRANDO: Como melhorar? RELEMBRANDO: A fórmula do método 3/8 de Simpson é Como melhorar? RELEMBRANDO: A fórmula do método 3/8 de Simpson é 3 minutos pra vocês modificarem o método 1/3 de Simpson para tornar-se o método 3/8 de Simpson. Go, go, go!!! 3/8 Simpson 3/8 Simpson Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; i = 0; Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; i = 0; for k = 0:n-1 Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (3/8)*((b-a)/3)*(f(a)+3*f((2*a+b)/3)+3*f((a+2*b)/3)+f(b)); Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (3/8)*((b-a)/3)*(f(a)+3*f((2*a+b)/3)+3*f((a+2*b)/3)+f(b)); end Scinotes 3/8 Simpson function i=simpson38(a0, bf, n) n=n/3; h=(bf-a0)/n; i = 0; for k = 0:n-1 a = a0+h*k; b = a0+h*(k+1); i = i + (3/8)*((b-a)/3)*(f(a)+3*f((2*a+b)/3)+3*f((a+2*b)/3)+f(b)); end endfunction Scinotes 3/8 Simpson Scilab 5.5.1 Console 3/8 Simpson -->simpson38(0,3,3) Scilab 5.5.1 Console 3/8 Simpson -->simpson38(0,3,3) ans = 1.3687612 Scilab 5.5.1 Console 3/8 Simpson -->simpson38(0,3,3) ans = 1.3687612 Scilab 5.5.1 Console 3/8 Simpson -->simpson38(0,3,3) ans = 1.3687612 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos. 3/8 Simpson -->simpson38(0,3,6) ans = 1.2767747 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos. 3/8 Simpson -->simpson38(0,3,9) ans = 1.2723078 Scilab 5.5.1 Console Os pontos extra são os azuis, e os vermelhos são os pontos dos N subintervalos.
Compartilhar