Métodos Numéricos e Série de Taylor: Aproximação de Funções

1 In33333 PONTIFÍCIA UNIVERSIDADE CATÓLICA DE CAMPINAS Engenharias CEATEC Introdução à Simulação Numérica Professores Alexandre Bia Cintia Denise e Vinícius Avaliações P1 ATIVIDADE P2 PROJETO Recuperação 1904 0604 2704 3105 1105 0806 1406 Unidade 6 MÉTODOS NUMÉRICOS DE RESOLUÇÃO DE EQUAÇÕES DIFERENCIAIS Exercícios para aula SÉRIES DE TAYLOR Considere uma função 𝑓𝑥 Em várias situações podemos utilizar uma aproximação dessa função por outras funções desde que a função satisfaça algumas condições O propósito dessa aproximação é utilizar funções cujas operações são mais simples como por exemplo funções polinomiais A aproximação de funções pode ser feita por meio da utilização de séries que são somas infinitas Uma das ferramentas matemáticas disponíveis para isso é a Série de Taylor que é uma série de potências em 𝑥 soma infinita de termos que envolvem potências da variável 𝑥 A série de Taylor analisa o comportamento da função 𝑓𝑥 quando 𝑥 está próximo de um valor fixo 𝑥0 Na prática quando a função 𝑓𝑥 admite uma representação por séries de potências o valor de 𝑓𝑥0 pode ser aproximado a partir da utilização de alguns termos da série de Taylor calculados em 𝑥0 Isto é o valor de 𝑓𝑥0 será aproximado pelo valor de um polinômio 𝑝𝑥 em 𝑥0 ou seja 𝑓𝑥0 𝑝𝑥0 Obviamente quando a função 𝑓𝑥 não é polinomial essa aproximação produz um erro que depende do grau do polinômio utilizado e pode ser estimado DEFINIÇÃO Se 𝑓𝑥 admite uma representação em série de potências do tipo 𝑓𝑥 𝑎𝑛𝑥 𝑥0𝑛 então 𝑓𝑘 existe para todo k ℕ e 𝑎𝑛 𝑓𝑛𝑥0 𝑛 Assim 𝑓𝑥 𝑓𝑥0 𝑓𝑥0𝑥 𝑥0 𝑓𝑥0 2 𝑥 𝑥02 𝑓𝑛𝑥0 𝑛 𝑥 𝑥0𝑛 A série acima é denominada SÉRIE DE TAYLOR A análise da série de Taylor fornece alguns resultados importantes Seja 𝑎𝑛 𝑥 𝑥0𝑛 𝑛0 uma série de Taylor Então OU i A série converge apenas para 𝑥 𝑥0 OU ii A série é convergente para qualquer valor de 𝑥0 OU iii A série converge num intervalo 𝑥 𝑥0 𝑟 O raio 𝑟 0 é chamado raio de convergência O polinômio 𝑓𝑘𝑥0 𝑘 𝑥 𝑥0𝑘 𝑛 𝑘0 é chamado Polinômio de Taylor de grau n Se a função 𝑓𝑥 possui 𝑛 1 derivadas num intervalo 𝐼 que contém 𝑥0 e se 𝑀 for uma cota superior para 𝑓𝑛1𝑥 em 𝐼 então 𝑅𝑛𝑥 𝑓𝑛1𝑥 𝑛 1 𝑥 𝑥0𝑛1 𝑀 𝑛 1 𝑥 𝑥0𝑛1 para 𝑥 𝐼 Exemplo 1 Encontre o polinômio de Taylor de grau 2 para as funções abaixo e use o polinômio para aproximar 2 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 o valor da função nos valores de 𝑐 indicados Compare com o valor exato a 𝑓𝑥 𝑠𝑒𝑛 𝑥 𝑥0 𝜋 3 𝑐 𝜋4 b 𝑓𝑥 𝑒𝑥 𝑥0 1 𝑐 05 Exemplo 2 Encontre a série de Taylor centrada em 𝑎 0 Considere que o intervalo de convergência da série de Taylor é Trace 𝑓 e seus primeiros 5 polinômios de Taylor na mesma tela O que você observa sobre a relação entre esses polinômios e 𝑓 a 𝑓𝑥 cos 𝑥 b 𝑓𝑥 𝑒𝑥 Solução 2a Arranjamos nossos cálculos em duas colunos como a seguir Como as derivadas se repetem em um ciclo de quatro podemos escrever a série de Taylor da seguinte forma 𝑓0 𝑓0 1 𝑥 𝑓0 2 𝑥2 𝑓0 3 1 0 1 𝑥 𝑥2 2 0𝑥3 3 𝑥4 4 1 𝑥2 2 𝑥4 4 𝑥6 6 cos𝑥 1𝑛𝑥2𝑛 2𝑛 𝑛0 para todo 𝑥 Vamos plotar os cinco primeiros polinômios de Taylor em cor azul e a função em cor vermelha na mesma janela usando o scilab Algoritmo para plotar os n primeiros polinômios da série de taylor e f variável auxiliar s0 domínio xpi001pi polinômios for i15 a 1i1factorial2i2x2i2 ssa plotxs end funcao cosseno spi001pi for i1lengths ficossi end plotsfr Exemplo 3 Use a série de Taylor de 𝑒𝑥 para calcular 𝑒02 Com precisão de 5 casas decimais sabendo que a série de Taylor de 𝑒𝑥 converge em Solução Primeiramente precisamos estimar quantos termos da série de Taylor precisaremos para conseguir calcular 𝑒02 pela série de Taylor com uma precisão de 104 Neste sentido vamos recorrer a fórmula do resto de Taylor 𝑅𝑛𝑥 𝑓𝑛1𝑥 𝑛 1 𝑥 0𝑛1 000001 𝑒𝑥 𝑛1 𝑥𝑛1 Escolhendo 𝑥 0 e atribuindo 𝑥 02 temos que 000001 𝑒0 𝑛1 02𝑛1 3 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Vamos calcular a contribuição do lado direito da equação usando um cálculo iterativo no scilab para estimar o valor de n que atenderia a equação for i110 i exp002i1factoriali1 end ans 1 002 2 00013333 3 00000667 4 00000027 5 8889D08 6 2540D09 7 6349D11 8 1411D12 9 2822D14 10 5131D16 Temos que precisaremos de 4 termos da série de Taylor para conseguir uma precisão de 5 casas decimais Então como 𝑒𝑥 𝑥𝑛 𝑛 𝑛0 temos que 𝑒𝑥 1 𝑥 𝑥2 2 𝑥3 3 Avaliando em 𝑥 02 nós temos que e02 1 02 022 2 023 3 e02 081866 APROXIMAÇÃO POR DIFERENÇAS FINITAS Aproximação por diferença centrada Se tivermos 3 pontos em uma malha 𝒙𝟎 𝒉 𝒙𝟎 𝒙𝟎 𝒉 Se você quiser calcular a derivada de uma função 𝑦 𝑓𝑥 ao longo do eixo x nós precisaremos desenvolver uma fórmula de diferenças entre os valores da função nestes três pontos que resultará em uma aproximação para a derivada da função no ponto 𝑥 Neste sentido iremos fazer uma aproximação de Taylor de 𝑓 em 𝑥 𝒙𝟎 de 3a ordem é dada por 𝑦𝑥 𝑦𝒙𝟎 𝑦𝒙𝟎𝑥 𝒙𝟎 𝑦𝒙𝟎𝑥 𝒙𝟎𝟐 2 𝑦𝒙𝟎𝑥 𝒙𝟎3 3 Avaliando a fórmula de Taylor em 𝑥 𝒙𝟎 𝒉 e em 𝑥 𝒙𝟎 𝒉 temos que 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉𝑦𝒙𝟎 𝒉2𝑦𝒙𝟎 2 𝒉3𝑦𝒙𝟎 3 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉𝑦𝒙𝟎 𝒉𝟐𝑦𝒙𝟎 2 𝒉3𝑦𝒙𝟎 3 Subtraindo as equações nós temos que 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉 𝟐𝒉𝑦𝒙𝟎 2𝒉3𝑦𝒙𝟎 6 Isolando 𝑦𝒙𝟎 temos que 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉 2𝒉3𝑦𝒙𝟎 6 𝟐𝒉𝑦𝒙𝟎 𝟐𝒉𝑦𝒙𝟎 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉 2𝒉3𝑦𝒙𝟎 6 𝑦𝒙𝟎 𝑦𝒙𝟎𝒉𝑦𝒙𝟎𝒉2𝒉3𝑦𝒙𝟎 6 2𝒉 Desmembrando apropriadamente o segundo membro temos que 𝑦𝒙𝟎 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉 2ℎ 𝒉2 𝑦𝒙𝟎 6 Dizemos que o termo 𝒉2 𝑦𝒙𝟎 6 é um termo de erro Ou seja desta forma obtemos uma fórmula discreta para a aproximação de 𝑦𝒙𝟎 usando os valores de 𝑦 nos pontos discretos das extremidades da malha 𝑦𝒙𝟎 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉 2ℎ 𝑂𝒉𝟐 Para obter a fórmula discreta de diferenças finitas de 𝑦𝒙𝟎 precisamos fazer um procedimento análogo ao anterior entretanto considerando uma expansão de Taylor de 4a ordem 4 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉𝑦𝒙𝟎 𝒉2𝑦𝒙𝟎 2 𝒉3𝑦𝒙𝟎 3 𝒉4𝑦𝒙𝟎 4 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉𝑦𝒙𝟎 𝒉2𝑦𝒙𝟎 2 𝒉3𝑦𝒙𝟎 3 𝒉4𝑦𝒙𝟎 4 Somando as duas equações acima temos que 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉 2𝑦𝒙𝟎 𝒉2𝑦𝒙𝟎 𝒉4 12 𝑦𝒙𝟎 Ou seja 𝒉2𝑦𝒙𝟎 𝑦𝒙𝟎 𝒉 𝑦𝒙𝟎 𝒉 2𝑦𝒙𝟎 𝒉4 12 𝑦𝒙𝟎 𝑦𝒙𝟎 𝑦𝒙𝟎 𝒉 2𝑦𝒙𝟎 𝑦𝒙𝟎 𝒉 𝒉2 𝑂ℎ2 Estas são as chamadas fórmulas de diferenças finitas centradas Elas serão especialmente úteis nos métodos de resolução numérica de equações diferenciais que iremos desenvolver Exemplo Aproxime o valor da 1ª e da 2ª derivada da função 𝑦 𝑒𝑥2 em x0 2 usando diferenças finitas com a h 01 b h 001 PROBLEMAS DE VALORES DE CONTORNO Exemplo Determine numericamente a altura atingida por uma pedra sabendose que após 5 segundos ela atinge a altura de 50 metros Solução Sabemos que a equação diferencial que descreve o movimento de queda livre é dada por 𝑑2𝑦 𝑑𝑡2 𝑔 Temos os seguintes dados 𝑔 98 𝑚 𝑠2 𝑦0 0 𝑚 𝑦5 50 𝑚 Vamos discretizar o domínio da variável t considerando 𝑛 1 pontos Ou seja se 𝑛 10 Então ℎ 5 0 10 ℎ 05 Então 𝑦0 0 𝑦𝑛 50 Discretizando a derivada segunda 𝑑2𝑦 𝑑𝑡2 a partir da fórmula de diferenças finitas centrada temos que 𝑑2𝑦 𝑑𝑡2 𝑦𝑖1 2𝑦𝑖 𝑦𝑖1 ℎ2 Ou seja substituindo esta contribuição na equação diferencial nós temos que 𝑦𝑖1 2𝑦𝑖 𝑦𝑖1 ℎ2 𝑔 𝑦𝑖1 2𝑦𝑖 𝑦𝑖1 𝑔ℎ2 𝑖 12 10 Usando a notação matricial temos que Observe que o problema de se resolver numericamente a equação diferencial se reduziu a um problema de se resolver um sistema linear com matriz dos coeficientes tridiagonal Abaixo apresentamos o código no scilab para construir a matriz tridiagonal do sistema linear em questão Algoritmo para criar a matriz tridiagonal n 11 a 2 b 1 c 1 A diagaones1n diagbones1n11 diagcones1n11 v1 1 zeros110 v2 zeros110 1 A1v1 A11v2 A A saída do programa será 5 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Vejamos o algoritmo para construir o vetor de termos independentes 𝒃 Algoritmo para construir o vetor de termos independentes g98 h05 b110 for i210 bi1 gh2 end b11150 b A saída do programa será Finalmente resolvemos o sistema linear usando a noção de inversão de matrizes 𝐴 𝒚 𝒃 Multiplica pela inversa de 𝐴 nos dois membros da equação acima 𝐴1 𝐴 𝒚 𝐴1 𝒃 𝐼 𝒚 𝐴1 𝒃 𝒚 𝐴1 𝒃 Ou seja a solução numérica da equação diferencial será a função 𝒚 𝑓𝒙 dada de forma discreta Fazendo o cálculo no scilab temos que Finalmente plotamos o gráfico de 𝒚 𝑓𝒙 nos pontos discretos x0055 xx yinvAb plotxy plotxyo EQUAÇÕES DIFERENCIAIS PARCIAIS Uma equação diferencial parcial EDP é uma equação envolvendo uma ou mais derivadas parciais de uma função de algumas variáveis que é desconhecida As equações diferenciais parciais que serão objeto de estudo estão associadas a três tipos distintos de fenômenos físicos processos de difusão processos oscilatórios e processos independentes do tempo ou estacionários Essas equações são portanto de importância fundamental em muitos ramos da física e da engenharia 6 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 EQUAÇÃO DO CALOR Inicialmente vamos analisar a seguinte equação diferencial parcial denominada equação de Poisson 2𝑢 𝑥2 𝑥 𝑦 2𝑢 𝑦2 𝑥 𝑦 𝑓𝑥 𝑦 em 𝑅 𝑥 𝑦𝑎 𝑥 𝑏 𝑐 𝑦 𝑑 com 𝑢𝑥 𝑦 𝑔𝑥 𝑦 para 𝑥 𝑦 𝑆 em que 𝑆 denota a fronteira de 𝑅 Se 𝑓 e 𝑔 forem contínuas em seu domínio então existe uma única solução para esta equação Escolha de uma malha O método utilizado é uma adaptação do método das diferenças finitas para problemas de contorno lineares justamente para a discretização do domínio O primeiro passo é escolher inteiros n e m e definir tamanhos de passo ℎ 𝑏𝑎 𝑛 e 𝑘 𝑑𝑐 𝑚 Divida o intervalo 𝑎 𝑏 em 𝑛 partes iguais de largura ℎ e intervalo 𝑐 𝑑 em 𝑚 partes iguais de largura 𝑘 conforme a figura abaixo Ponha uma grade no retângulo 𝑅 desenhando linhas verticais e horizontais através dos pontos com coordenadas 𝑥𝑖 𝑦𝑗 em que 𝑥𝑖 𝑎 𝑖ℎ para 𝑖 01 𝑛 𝑦𝑗 𝑏 𝑗𝑘 para 𝑗 01 𝑚 As retas 𝑥 𝑥𝑖 e 𝑦 𝑦𝑗 são as retas da grade e suas intersecções são os pontos da malha da grade Para cada ponto da malha no interior da grade 𝑥𝑖 𝑦𝑗 para 𝑖 12 𝑛 1 e 𝑗 12 𝑚 1 Usamos as fórmulas de diferenças finitas centradas que podem ser deduzidas a partir da expansão em série de Taylor a exemplo de como foi demonstrado no início da unidade 2𝑢 𝑥2 𝑥𝑖 𝑦𝑗 𝑢𝑥𝑖1 𝑦𝑗 2𝑢𝑥𝑖 𝑦𝑗 𝑢𝑥𝑖1 𝑦𝑗 ℎ2 2𝑢 𝑦2 𝑥𝑖 𝑦𝑗 𝑢𝑥𝑖 𝑦𝑗1 2𝑢𝑥𝑖 𝑦𝑗 𝑢𝑥𝑖 𝑦𝑗1 𝑘2 Substituindo estas contribuições na equação de Poisson temos que 2𝑢 𝑥2 𝑥 𝑦 2𝑢 𝑦2 𝑥 𝑦 𝑓𝑥 𝑦 𝑢𝑥𝑖1𝑦𝑗2𝑢𝑥𝑖𝑦𝑗𝑢𝑥𝑖1𝑦𝑗 ℎ2 𝑢𝑥𝑖𝑦𝑗12𝑢𝑥𝑖𝑦𝑗𝑢𝑥𝑖𝑦𝑗1 𝑘2 𝑓𝑥𝑖 𝑦𝑗 Definimos 𝑤𝑖𝑗 𝑢𝑥𝑖 𝑦𝑗 e então temos que 𝑤𝑖1𝑗 2𝑤𝑖𝑗 𝑤𝑖1𝑗 ℎ2 𝑤𝑖𝑗1 2𝑤𝑖𝑗 𝑤𝑖𝑗1 𝑘2 𝑓𝑥𝑖 𝑦𝑗 Multiplicando a equação acima por ℎ2 nós temos que 𝑤𝑖1𝑗 2𝑤𝑖𝑗 𝑤𝑖1𝑗 2𝑤𝑖𝑗 ℎ2 𝑘2 𝑤𝑖𝑗1 ℎ2 𝑘2 𝑤𝑖𝑗1 ℎ2 𝑘2 ℎ2 𝑓𝑥𝑖 𝑦𝑗 Reorganizando os termos temos que 2 ℎ 𝑘 2 1 𝑤𝑖𝑗 𝑤𝑖1𝑗 𝑤𝑖1𝑗 ℎ2 𝑘2 𝑤𝑖𝑗1 𝑤𝑖𝑗1 ℎ2 𝑓𝑥𝑖 𝑦𝑗 I Quanto as condições de contorno temos que 𝑤0𝑗 𝑔𝑥0 𝑦𝑗 e 𝑤𝑛𝑗 𝑔𝑥𝑛 𝑦𝑗 𝑤𝑖0 𝑔𝑥𝑖 𝑦0 e 𝑤𝑖𝑚 𝑔𝑥𝑖 𝑦𝑚 A Equação I envolve aproximações de 𝑢𝑥 𝑦 nos pontos 𝑥𝑖1 𝑦𝑗 𝑥𝑖 𝑦𝑗 𝑥𝑖1 𝑦𝑗 𝑥𝑖 𝑦𝑗1 e 𝑥𝑖 𝑦𝑗1 A reprodução da porção da grade onde estes pontos estão localizados mostra que cada equação envolve aproximações em uma região em formato de estrela ao redor de 𝑥𝑖 𝑦𝑗 Conforme exibe a figura abaixo 7 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Usamos a informação das condições de contorno sempre que apropriado no sistema dado pela Equação I isto é em todos os pontos 𝑥𝑖 𝑦𝑗 adjacentes a um ponto da fronteira da malha Isso produz um sistema linear 𝑛 1𝑚 1 𝑛 1𝑚 1 com as incógnitas sendo as aproximações 𝑤𝑖𝑗 para 𝑢𝑥𝑖 𝑦𝑗 nos pontos interiores da malha Para expressar os cálculos com matriz de forma mais eficiente vamos introduzir uma enumeração dos pontos interiores da malha 𝑃𝑙 𝑥𝑖 𝑦𝑗 e 𝑤𝑙 𝑤𝑖𝑗 Em que 𝑙 𝑖 𝑚 1 𝑗𝑛 1 para cada 𝑖 12 𝑛 1 e 𝑗 12 𝑚 1 Isto enumera os pontos da malha consecutivamente da esquerda para a direita e de cima para baixo Enumerar os pontos desta maneira assegura que o sistema necessário para determinar 𝑤𝑖𝑗 diga respeito a uma matriz de banda com largura máxima de 2𝑛 1 Exemplo Determine a distribuição de calor estacionária em uma placa de metal fina com dimensões 05𝑚 por 05𝑚 usando 𝑛 𝑚 4 Dois lados adjacentes são mantidos a 0 𝐶 1 𝑜 e o calor nos outros aumenta linearmente de 0 𝐶 1 𝑜 em um canto a 100 𝐶 1 𝑜 onde os lados se encontram Solução Se colocarmos os lados com condições de contorno nulas ao longo dos eixos x e y o problema será expressado essencialmente como 2𝑢 𝑥2 𝑥 𝑦 2𝑢 𝑦2 𝑥 𝑦 0 para cada 𝑥 𝑦 em 𝑅 𝑥 𝑦0 𝑥 05 0 𝑦 05 com as condições de contorno 𝑢0 𝑦 0 𝑢𝑥 0 0 𝑢𝑥 05 200𝑥 e 𝑢05 𝑦 200𝑦 O problema tem a grade abaixo Os parâmetros ℎ e 𝑘 são calculados assim ℎ 𝑏 𝑎 𝑛 05 0 4 ℎ 0125 𝑘 𝑑 𝑐 𝑘 05 0 4 𝑘 0125 Substituindo os dados na Equação I temos que 2 ℎ 𝑘 2 1 𝑤𝑖𝑗 𝑤𝑖1𝑗 𝑤𝑖1𝑗 ℎ2 𝑘2 𝑤𝑖𝑗1 𝑤𝑖𝑗1 ℎ2 𝑓𝑥𝑖 𝑦𝑗 2 0125 0125 2 1 𝑤𝑖𝑗 𝑤𝑖1𝑗 𝑤𝑖1𝑗 01252 01252 𝑤𝑖𝑗1 𝑤𝑖𝑗1 01252 0 Ou seja 4𝑤𝑖𝑗 𝑤𝑖1𝑗 𝑤𝑖1𝑗 𝑤𝑖𝑗1 𝑤𝑖𝑗1 0 para cada 𝑖 123 e 𝑗 123 Se expressarmos as equações em termos dos pontos interiores da grade lembrando que 𝑤𝑖 𝑢𝑃𝑖 e 𝑤𝑙 𝑤𝑖𝑗 teremos que Por exemplo de acordo com a malha o ponto 𝑃1está localizado pelos índices 𝑖 1 e 𝑗 3 Com efeito 𝑙 𝑖 𝑚 1 𝑗𝑛 1 𝑙 1 4 1 34 1 𝑙 1 0 3 𝑙 1 Neste ponto 𝑃1teremos que 8 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 4𝑤𝑖𝑗 𝑤𝑖1𝑗 𝑤𝑖1𝑗 𝑤𝑖𝑗1 𝑤𝑖𝑗1 0 4𝑤13 𝑤113 𝑤113 𝑤131 𝑤131 0 4𝑤13 𝑤23 𝑤03 𝑤14 𝑤12 0 Usando o fato que 𝑤𝑙 𝑤𝑖𝑗 teremos que 4𝑤1 𝑤2 𝑤03 𝑤14 𝑤4 0 Reescrevendo a equação em termos das variáveis 𝑤𝑙 temos que 𝑃1 4𝑤1 𝑤2 𝑤4 𝑤03 𝑤14 Trabalhando similarmente para os demais pontos temos que em que os lados direitos das equações são obtidos a partir das condições de contorno De fato as condições de contorno implicam que 𝑤10 𝑤20 𝑤30 𝑤01 𝑤02 𝑤03 0 𝑤14 𝑤41 25 𝑤24 𝑤42 50 e 𝑤34 𝑤43 75 Neste sentido o sistema linear associado a este problema tem a forma Observe que a matriz dos coeficientes do sistema linear resultou em uma matriz banda esparsa simétrica e tridiagonal por blocos onde Vamos resolver este problema no scilab Primeiramente vamos construir o bloco tridiagonal depois o bloco da oposta da identidade e finalmente construir nossa matriz dos coeficientes definida por blocos N3 a4 b1 c1 A diagaonesN1diagbones1N 11diagcones1N11 Beye33 Czeros33 M A B CB A BC B A Agora vamos resolver o sistema linear 𝑀 𝑥 𝑏 M A B CB A BC B A b 25 50 150 0 0 50 0 0 25 xinvMb Saída do algoritmo Plotando a superfície da solução x 0012505 9 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 y 0012505 Zzeros15 0 625 125 1875 250 125 25 375 500 1875 375 5625 75 0 25 50 75 100 meshxyZ EQUAÇÃO DA CORDA Modelamento da Equação da Corda Os problemas de vibrações em cordas são muito estudados pelos físicos e pelos engenheiros em aplicações críticas Uma representação deste problema é mostrada na figura a seguir A EDP que modela o comportamento dinâmico de uma corda é dada por 2𝑦 𝑥2 𝜇 𝑇 2𝑦 𝑡2 𝑘 𝑦 𝑡 1 Neste caso 𝑦𝑥 𝑡 é a posição vertical de um ponto 𝑥 da corda no instante 𝑡 𝜇 é a densidade mássica por comprimento da corda 𝑇 é a tração na corda e 𝑘 é um coeficiente de dissipação de energia Aplicação do Método das Diferenças Finitas Para aplicar o método das diferenças finitas é necessário primeiramente definir a notação utilizada Tanto no domínio do espaço 𝑥 quanto no domínio do tempo 𝑡 propõese uma discretização das variáveis independentes A variável independente 𝑥 será discretizada de forma a ser dividida em 𝑁𝑥 segmentos de igual tamanho ℎ𝑥 𝐿𝑁𝑥 Assim o cálculo é feito no espaço sobre 𝑁𝑥 1 variáveis sendo que para 𝑥 0 o índice em 𝑥 da variável é 0 e em 𝑥 𝐿 o índice em 𝑥 da variável é 𝑁𝑥 A variável independente 𝑡 é discretizada de maneira análoga em 𝑁𝑡 segmentos de tamanho ℎ𝑡 𝑇𝑓 𝑁𝑡 sendo 𝑇𝑓 o tempo final da análise Assim para o tempo são 𝑁𝑡 1 segmentos iniciando em 0 e terminando em 𝑁𝑡 Assim definese a notação para a função 𝑦𝑥0 𝑖ℎ𝑥 𝑡0 𝑗ℎ𝑡 𝑦𝑖𝑗 Além disso se a própria variável independente for mencionada no subíndice esta notação representa um conjunto de equações para qualquer valor daquela variável Exemplos para condição inicial podese escrever 𝑦𝑥0 0 o que representa que 𝑦𝑥 𝑡 0 0 para qualquer valor de 𝑥 Para condição de contorno podese escrever 𝑦𝐿𝑡 0 o que significa que 𝑦𝑥 𝐿 𝑡 0 para qualquer valor de 𝑡 Com estas considerações já é possível calcular as condições de contorno do problema Estas condições para o problema da corda são simples pois as extremidades da corda não podem se mover na vertical Assim podese escrever 𝑦0𝑡 0 𝑦𝑁𝑥𝑡 0 2 Já para as condições iniciais partese do pressuposto que a posição e a velocidade de todos os pontos da corda são conhecidas no tempo 𝑡 0 Assim podese escrever 𝑦𝑥0 𝑦0𝑥 𝑑𝑦𝑥0 𝑑𝑡 𝑦0𝑥 3 No caso da segunda condição inicial da equação 3 podese aplicar a ela a aproximação da 10 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 primeira derivada da função 𝑦𝑥0 com relação a 𝑡 como sendo 𝑑𝑦𝑥0 𝑑𝑡 𝑦𝑥1𝑦𝑥1 2ℎ𝑡 𝑦0𝑥 4 O que resulta na relação 𝑦 𝑥1 𝑦𝑥1 2ℎ𝑡𝑦0𝑥 5 Pelo fato de se estar aplicando o método das diferenças finitas centrado no tempo é necessário também criar uma condição de encerramento da integração numérica Esta condição é realizada através da aproximação do último passo de integração no tempo utilizando a aproximação de primeira ordem da integral de Newton Assim tem se 𝑦𝑥𝑁𝑡 𝑦𝑥𝑁𝑡1 ℎ𝑡 𝑑𝑦𝑥𝑁𝑡1 𝑑𝑥 6 Novamente é possível aproximar a derivada de 𝑦𝑥𝑁𝑡1 pelo método das diferenças finitas Assim se obtém 𝑦𝑥𝑁𝑡 𝑦𝑥𝑁𝑡1 ℎ𝑡 𝑦𝑥𝑁𝑡 𝑦𝑥𝑁𝑡2 2ℎ𝑡 7 Aplicando as simplificações cabíveis e reorganizando se obtém a equação 𝑦𝑥𝑁𝑡 2𝑦𝑥𝑁𝑡1 𝑦𝑥𝑁𝑡2 0 8 A equação obtida em 8 permite realizar a finalização da montagem da matriz do método das diferenças finitas sem a necessidade de uma nova condição de contorno ou de tempo Assim as condições de contorno iniciais e de finalização são reunidas na equação 9 𝑦0𝑡 0 𝑦𝑁𝑥𝑡 0 𝑦𝑥0 𝑦0𝑥 𝑦𝑥1 𝑦𝑥1 2ℎ𝑡𝑦0𝑥 𝑦𝑥𝑁𝑡 2𝑦𝑥𝑁𝑡1 𝑦𝑥𝑁𝑡2 0 9 Além de todas estas condições ainda é necessário obter a regra geral do sistema que vem da aplicação do método das diferenças finitas à equação 1 Realizando esta aplicação se obtém 𝑦𝑖1𝑗 𝑦𝑖1𝑗 2𝑦𝑖𝑗 ℎ𝑥2 𝜇 𝑇 𝑦𝑖𝑗1 𝑦𝑖𝑗1 2𝑦𝑖𝑗 ℎ𝑡 2 𝑘 𝑦𝑖𝑗1 𝑦𝑖𝑗1 2ℎ𝑡 0 10 Notase que a equação 10 possui índices de 𝑖 1 a 𝑖 1 no espaço e de 𝑗 1 a 𝑗 1 no tempo As condições de contorno permitem calcular os elementos 0 e 𝑁𝑥 no espaço e as condições temporais permitem o cálculo dos índices 1 0 e 𝑁𝑡 no tempo Assim a regra geral deve ser aplicada para 1 𝑖 𝑁𝑥 1 e 1 𝑗 𝑁𝑡 1 Utilizando as relações das equações 9 e 10 devese agora escrever o problema na forma 𝐴𝑤 𝑞 e inverter a matriz 𝐴 para obter a solução desejada O PROBLEMA DA VIGA 1 Modelamento Matemático da Viga O problema de flexão de vigas é um problema amplamente estudado na engenharia Tratase de estudar um problema no qual se calcula os esforços internos e deflexões de vigas resultantes de suas propriedades geometria e dos carregamentos aplicados nela A viga é um elemento estrutural delgado caracterizado por suportar esforços perpendiculares ao seu eixo Outros elementos estruturais importantes são os eixos que suportam torção axial e barras que suportam esforços axiais 11 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Há dois casos de vinculação condições de contorno muito estudados na engenharia que são as vigas biapoiadas e as vigas em balanço conforme ilustração a seguir As vigas como mencionado anteriormente são sujeitas a cargas perpendiculares a seus eixos conforme a figura a seguir O eixo 𝑥 que é a variável independente de análise da viga no espaço possui sua origem no início da viga normalmente no lado esquerdo O eixo 𝑦 por convenção aponta para baixo pois normalmente é neste sentido que são aplicadas as cargas em uma viga e sua deflexão neste sentido é comumente positiva O carregamento aplicado à viga é uma função descrita em 𝑥 normalmente denominada 𝑞𝑥 11 Caso estático O caso estático é o caso em que todas as oscilações temporais na resposta da viga já cessaram e não é mais necessário considerar variações temporais Este caso é importante na construção da EDP que descreve a viga pois servirá de base para o caso dinâmico com as oscilações temporais O esforço cortante é um esforço interno cisalhante que provoca tensões perpendiculares ao eixo da viga denominado por 𝑉𝑥 Este esforço é dado pela expressão 𝑑𝑉𝑥 𝑑𝑥 𝑞𝑥 1 O momento fletor 𝑀𝑥 também é um esforço interno da viga mas este causa tensões paralelas ao eixo da viga Ele é dado por 𝑑𝑀𝑥 𝑑𝑥 𝑉𝑥 2 O ângulo de deflexão da viga 𝜃𝑥 é um resultado direto dos esforços internos da viga de seu material representado pelo módulo de Young 𝐸 e de sua geometria representada pelo momento de inércia de área 𝐼𝑧𝑧 Ele é calculado através da seguinte equação 𝐸𝐼𝑧𝑧 𝑑𝜃𝑥 𝑑𝑥 𝑀𝑥 3 Já a deflexão ou flecha da viga 𝑤𝑥 é calculada com base no ângulo de deflexão provocado pelos esforços internos através da seguinte equação 𝑑𝑤𝑥 𝑑𝑥 𝜃𝑥 4 Unindo os resultados encontrados nas equações de 1 a 4 é possível construir a EDO que relaciona a flecha da viga e o carregamento a ela aplicado no caso estático 𝐸𝐼𝑧𝑧 𝑑4𝑤𝑥 𝑑𝑥4 𝑞𝑥 5 Esta equação é uma EDO pois se trata de um caso estático no qual as oscilações temporais são desconsideradas e o carregamento não varia no tempo Contudo no caso de variações temporais esta equação se transformará em uma EDP cujas variáveis independentes serão o comprimento 𝑥 e o tempo 𝑡 Este caso será tratado mais adiante Antes de avançar para o caso dinâmico é importante a partir do caso estático definir as condições de contorno do problema para os casos de viga biapoiada e em balanço Para isso é importante compreender não apenas a EDO mas também as equações de 1 a 4 que regem o problema 12 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Como convenção adotase que a extremidade esquerda da viga se encontra em 𝑥 0 e a extremidade direita em 𝑥 𝐿 sendo 𝐿 o comprimento da viga No caso da viga biapoiada o deslocamento vertical em ambas as extremidades é nulo Assim o primeiro par de condições de contorno é dado por 𝑤𝑥 0 0 e 𝑤𝑥 𝐿 0 Os vínculos no caso da viga biapoiada são capazes de suportar os esforços verticais mas não os momentos fletores Assim para que a viga esteja em equilíbrio o segundo par de condições de contorno é 𝑀𝑥 0 0 e 𝑀𝑥 𝐿 0 Já no caso da viga em balanço considerando o engaste na extremidade esquerda da viga a flecha e o ângulo de deflexão nesta extremidade serão nulos Assim temse o primeiro par de condições de contorno que é 𝑤𝑥 0 0 e 𝜃𝑥 0 0 Na extremidade esquerda não há nenhum vínculo para equilibrar os esforços internos assim o momento fletor e o esforço cortante nesta extremidade devem ser nulos Isto fornece o segundo par de condições de contorno 𝑀𝑥 𝐿 0 e 𝑉𝑥 𝐿 0 Estas condições de contorno vão se relacionar diretamente com a aplicação do método das diferenças finitas que se deseja aplicar para a solução do problema e são válidas também para o caso dinâmico no qual além das condições de contorno é necessário também obter condições iniciais para o problema 12 Caso dinâmico No caso dinâmico além de todas as considerações do caso estático considerase também as variações temporais da resposta da viga Isto significa que a carga aplicada pode variar com o tempo e que pelo fato de a viga possuir inércia e a capacidade de dissipar energia em deformação além de rigidez estática seu comportamento pode apresentar respostas transientes Respostas transientes são respostas que correspondem à adaptação do sistema à presença de um novo carregamento Na presença de dissipação de energia este transiente tende a zero e a resposta que se obtém após o tempo de resposta transiente é chamada de resposta permanente Para cargas invariantes no tempo a resposta permanente é exatamente a resposta do caso estático A EDP que representa o caso dinâmico é dada por 𝐸𝐼𝑧𝑧 4𝑤𝑥 𝑡 𝑥4 𝐴𝜌 2𝑤𝑥 𝑡 𝑡2 𝑘 𝑤𝑥 𝑡 𝑡 𝑞𝑥 𝑡 6 Neste caso 𝐴 é a área da seção transversal da viga considerada constante neste caso em 𝑥 𝜌 é a densidade volumétrica do material constante em toda a viga e 𝑘 é um coeficiente de dissipação de energia normalmente pequeno Para o caso dinâmico além das condições de contorno é necessário definir as condições iniciais do problema Como o grau de diferenciação no tempo é o grau 2 então são necessárias duas condições de contorno Uma delas é a posição inicial de todos os elementos da viga que pode ser descrita por 𝑤𝑥 𝑡 0 𝑤0𝑥 A segunda condição inicial é a velocidade dos pontos da viga que pode ser descrita por 𝑑𝑤𝑥 𝑡 0𝑑𝑡 𝑣0𝑥 1 Aplicação do Método das Diferenças Finitas Conforme visto nas seções anteriores a aproximação da primeira derivada de uma função genérica 𝑦𝑥 em função de 𝑥 em torno de um ponto genérico 𝑥0 é dada por 𝑑𝑦𝑥 𝑑𝑥 𝑦𝑥0 ℎ 𝑦𝑥0 ℎ 2ℎ 7 Esta aproximação é comumente referida como sendo a aproximação centrada por diferenças finitas da primeira derivada de 𝑦𝑥 A aproximação da segunda derivada desta função também é conhecida 𝑑2𝑦𝑥 𝑑𝑥2 𝑦𝑥0 ℎ 𝑦𝑥0 ℎ 2𝑦𝑥0 ℎ2 8 13 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 O leitor que acompanhou as deduções destas duas aproximações pode notar que através do mesmo método não é possível obter as aproximações das derivadas de ordens maiores Isto se deve ao fato de que nas aproximações de primeira e segunda ordem já se assumiu que os erros de ordem 2 e 3 são desprezíveis Assim as expansões em série de Taylor de ordem mais alta resultam apenas na igualdade trivial 0 0 Assim para se obter uma aproximação para a derivada de ordem 3 da função 𝑦𝑥 em torno de um ponto 𝑥0 é necessário derivar uma vez a aproximação de segunda ordem Assim 𝑑3𝑦𝑥 𝑑𝑥3 𝑑𝑦 𝑑𝑥 𝑥0 ℎ 𝑑𝑦 𝑑𝑥 𝑥0 ℎ 2 𝑑𝑦 𝑑𝑥 𝑥0 ℎ2 9 Utilizando a equação 7 para aproximar a primeira derivada de 𝑦𝑥 temse 𝑑3𝑦𝑥 𝑑𝑥3 𝑦𝑥0 2ℎ 𝑦𝑥0 2ℎ 𝑦𝑥0 𝑦𝑥0 2ℎ 2ℎ 2 𝑦𝑥0 ℎ 𝑦𝑥0 ℎ 2ℎ ℎ2 10 Portanto 𝑑3𝑦𝑥 𝑑𝑥3 𝑦𝑥0 2ℎ 2𝑦𝑥0 ℎ 2𝑦𝑥0 ℎ 𝑦𝑥0 2ℎ 2ℎ3 11 Para simplificar a notação adotase 𝑦𝑥0 𝑖ℎ 𝑦𝑖 Assim a equação 11 pode ser reescrita como 𝑑3𝑦𝑥 𝑑𝑥3 𝑦22𝑦12𝑦1𝑦2 2ℎ3 12 O mesmo procedimento pode ser feito para aproximar a quarta derivada da função 𝑦𝑥 gerando a seguinte equação 𝑑4𝑦𝑥 𝑑𝑥4 𝑑2𝑦 𝑑𝑥2𝑥0ℎ𝑑2𝑦 𝑑𝑥2𝑥0ℎ2𝑑2𝑦 𝑑𝑥2𝑥0 ℎ2 13 Então fazendo uso da equação 8 a equação 13 se transforma em 𝑑4𝑦𝑥 𝑑𝑥3 𝑦2 𝑦0 2𝑦1 𝑦0 𝑦2 2𝑦1 2𝑦1 𝑦1 2𝑦0 ℎ4 14 Ou seja 𝑑4𝑦𝑥 𝑑𝑥4 𝑦24𝑦16𝑦04𝑦1𝑦2 ℎ4 15 11 Caso Estático Viga Biapoiada Para a viga biapoiada as condições de contorno são 𝑤𝑥 0 0 𝑤𝑥 𝐿 0 𝑀𝑥 0 0 𝑀𝑥 𝐿 0 16 Para o caso de uma viga podese discretizar a função 𝑤𝑥 em 𝑁 1 pontos de forma que ℎ 𝐿𝑁 Assim o primeiro ponto é o ponto de índice 0 e fica em 𝑥 0 enquanto o último ponto fica em 𝑥 𝐿 e possui índice 𝑁 As duas primeiras condições de contorno da equação 16 são diretamente aplicáveis pois resultam em 𝑤0 0 e 𝑤𝑁 0 Já as duas últimas condições de contorno necessitam de atenção especial Combinando as equações 3 e 4 obtémse a relação 𝐸𝐼𝑧𝑧 𝑑2𝑤𝑥 𝑑𝑥2 𝑀𝑥 17 Aplicando a equação 8 na aproximação da segunda derivada de 𝑤𝑥 se obtém para a terceira condição de contorno 14 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 𝐸𝐼𝑧𝑧 𝑤1𝑤12𝑤0 ℎ2 𝑀𝑥 0 0 18 Neste caso 𝑤0 é dado pela condição de contorno anterior e 𝑤1 é uma variável do sistema a se calcular Assim devese obter uma equação para 𝑤1 𝑤1 2𝑤0 𝑤1 19 Já para a quarta condição de contorno se obtém 𝐸𝐼𝑧𝑧 𝑤𝑁1 𝑤𝑁1 2𝑤𝑁 ℎ2 𝑀𝑥 𝐿 0 20 Neste caso o 𝑤𝑁 é dado por uma condição de contorno anterior e o termo 𝑤𝑁1 é uma variável do sistema Assim devese encontrar uma equação para o termo 𝑤𝑁1 𝑤𝑁1 2𝑤𝑁 𝑤𝑁1 21 Assim as condições iniciais forneceram relações para os termos 𝑤1 𝑤0 𝑤𝑁 e 𝑤𝑁1 Para os termos 𝑤1 a 𝑤𝑁1 é possível utilizar a EDO do caso estático mostrada na equação 5 Utilizando esta equação e a aproximação da equação 15 é possível obter a relação 𝐸𝐼𝑧𝑧 𝑤𝑖2 4𝑤𝑖1 6𝑤𝑖 4𝑤𝑖1 𝑤𝑖2 ℎ4 𝑞𝑖 𝑝𝑎𝑟𝑎 1 𝑖 𝑁 1 22 Note que para 𝑖 1 o termo 𝑤 de menor índice acessado será o termo 𝑤1 e para 𝑖 𝑁 1 o termo 𝑤 de maior índice acessado será o termo 𝑤𝑁1 Ambos os termos já estão definidos pelas condições de contorno do sistema O termo 𝑞𝑖 é definido como 𝑞𝑥0 𝑖ℎ Como o carregamento da viga é conhecido basta calcular este termo para as várias posições desejadas Assim o sistema pode ser escrito na forma 𝐴𝑤 𝑞 sendo 𝑤 𝑤1 𝑤0 𝑤1 𝑤𝑁1 𝑤𝑁 𝑤𝑁1 𝑇 23 𝑞 0 0 𝑞1 𝑞𝑁1 0 0 𝑇 24 𝐴 1 2 1 0 0 0 0 0 0 1 0 0 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 0 0 1 0 0 0 0 0 0 1 2 1 25 Na equação 25 𝑝 𝐸𝐼𝑧𝑧ℎ4 12 Caso Estático Viga em Balanço Iniciase novamente por transformar as condições de contorno em equações que farão parte da solução do sistema As condições de contorno para este caso são 𝑤𝑥 0 0 𝜃𝑥 0 0 𝑀𝑥 𝐿 0 𝑉𝑥 𝐿 0 26 A primeira condição de contorno é a mais simples pois resulta diretamente em 𝑤0 0 A relação resultante da segunda condição de contorno pode ser calculada a partir da equação 4 e da equação 7 𝜃𝑥 0 𝑑𝑤𝑥0 𝑑𝑥 𝑤1𝑤1 2ℎ 0 27 De onde se extrai a relação de que 𝑤1 𝑤1 0 Para a terceira condição de contorno podese utilizar a equação 17 em conjunto com a equação 8 para obter 𝑀𝑥 𝐿 𝐸𝐼𝑧𝑧 𝑑2𝑤𝑥𝐿 𝑑𝑥2 𝐸𝐼𝑧𝑧 𝑤𝑁1𝑤𝑁12𝑤𝑁 ℎ2 0 28 De onde se conclui que 𝑤𝑁1 2𝑤𝑁 𝑤𝑁1 0 15 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Finalmente da última condição de contorno utilizase as equações 2 12 e 17 para obter 𝑉𝑥 𝐿 𝐸𝐼𝑧𝑧 𝑑3𝑤𝑥𝐿 𝑑𝑥3 𝐸𝐼𝑧𝑧 𝑤𝑁22𝑤𝑁12𝑤𝑁1𝑤𝑁2 2ℎ3 0 29 De onde se conclui que 𝑤𝑁2 2𝑤𝑁1 2𝑤𝑁1 𝑤𝑁2 0 Neste caso enquanto a primeira condição de contorno cria uma relação para 𝑤0 e a segunda uma relação para 𝑤1 a terceira cria uma relação para 𝑤𝑁1 e a quarta uma relação para 𝑤𝑁2 Assim neste caso específico a relação dada pela equação 22 deve ser usada no cálculo de 𝑤1 até 𝑤𝑁 e não 𝑤𝑁1 como no caso da viga biapoiada Desta forma o problema pode novamente ser escrito na forma 𝐴𝑤 𝑞 sendo 𝑤 𝑤1 𝑤0 𝑤1 𝑤𝑁 𝑤𝑁1 𝑤𝑁2 𝑇 30 𝑞 0 0 𝑞1 𝑞𝑁 0 0 𝑇 31 𝐴 1 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 𝑝 4𝑝 6𝑝 4𝑝 𝑝 0 0 0 0 1 2 1 0 0 0 0 1 2 0 2 1 3 2 13 Caso Dinâmico Viga Biapoiada Nos casos dinâmicos a solução será a de uma EDP e não a de uma EDO Por isso devese levar em consideração que a discretização das variáveis independentes será feita neste caso em dois domínios no comprimento e no tempo Assim no comprimento discretizase a variável 𝑥 em 𝑁𝑥 segmentos de igual comprimento ℎ𝑥 e no tempo se discretiza a variável 𝑡 em 𝑁𝑡 segmentos de igual comprimento ℎ𝑡 A variável 𝑤𝑥 𝑡 agora pode ser aproximada em torno de um ponto 𝑥0 𝑡0 de forma que para facilitar a escrita definese a apresentação 𝑤𝑥0 𝑖ℎ𝑥 𝑡0 𝑗ℎ𝑡 𝑤𝑖𝑗 Além disso quando o nome da variável independente for mencionado no subíndice então entendese um conjunto de equações para cada valor daquela variável Ex 𝑤𝑥0 0 representa que o valor de 𝑤𝑥 𝑡 0 0 para qualquer valor de 𝑥 Alternativamente 𝑤𝑝𝑡 0 representa que o valor de 𝑤𝑥 𝑝 𝑡 0 para qualquer valor de 𝑡 Novamente as condições de contorno deste problema são as mesmas da equação 16 repetida abaixo com a dependência do tempo 𝑤𝑥 0 𝑡 0 𝑤𝑥 𝐿 𝑡 0 𝑀𝑥 0 𝑡 0 𝑀𝑥 𝐿 𝑡 0 33 Assim para as duas primeiras condições de contorno temse 𝑤0𝑡 0 e 𝑤𝑁𝑥𝑡 0 Para a terceira temse 𝐸𝐼𝑧𝑧 𝑤1𝑡𝑤1𝑡2𝑤0𝑡 ℎ𝑥2 𝑀𝑥 0 0 34 De onde se conclui que 𝑤1𝑡 𝑤1𝑡 2𝑤0𝑡 0 Finalmente para a quarta condição de contorno temse 𝐸𝐼𝑧𝑧 𝑤𝑁𝑥1𝑡𝑤𝑁𝑥1𝑡2𝑤𝑁𝑥𝑡 ℎ𝑥2 𝑀𝑥 𝐿 0 35 Resultando em 𝑤𝑁𝑥1𝑡 𝑤𝑁𝑥1𝑡 2𝑤𝑁𝑥𝑡 0 Para as condições iniciais pressupõese o conhecimento dos deslocamentos e das velocidades iniciais em qualquer ponto da viga Assim para o deslocamento 𝑤𝑥0 𝑤0𝑥 e para a velocidade 16 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 𝑑𝑤𝑥0 𝑑𝑡 𝑤𝑥1𝑤𝑥1 2ℎ𝑡 𝑤 0𝑥 36 Assim 𝑤𝑥1 𝑤𝑥1 2ℎ𝑡𝑤 0𝑥 37 Para todas as outras variáveis 𝑤𝑖𝑗 não definidas pelas condições iniciais ou de contorno podese utilizar a equação 6 como a regra para a montagem da matriz para a simulação do sistema Aplicando o método das diferenças finitas sobre esta equação obtémse 𝐸𝐼𝑧𝑧 𝑤𝑖2𝑗 4𝑤𝑖1𝑗 6𝑤𝑖𝑗 4𝑤𝑖1𝑗 𝑤𝑖2𝑗 ℎ𝑥4 𝐴𝜌 𝑤𝑖𝑗1𝑤𝑖𝑗12𝑤𝑖𝑗 ℎ𝑡 2 𝑘 𝑤𝑖𝑗1𝑤𝑖𝑗1 2ℎ𝑡 𝑞𝑖𝑗 38 38 No espaço as condições de contorno são capazes de fornecer informações para a solução dos índices 1 0 𝑁𝑥 e 𝑁𝑥 1 Como a regra geral da equação 38 define que no espaço trabalhase com os coeficientes de 𝑖 2 a 𝑖 2 então no espaço a regra geral deve ser aplicada nos coeficientes de 1 a 𝑁𝑥 1 Já no tempo as condições iniciais fornecem subsídios para o cálculo dos índices 0 e 1 Como na regra geral equação 38 trabalhase entre os índices 𝑗 1 e 𝑗 1 então no tempo devese aplicar a regra para os índices de 1 a 𝑁𝑡 1 Assim o número de variáveis a se resolver é de 𝑁𝑡 2𝑁𝑥 3 Além de todas estas condições quando se trata da solução do sistema no tempo o método não propõe uma regra para o cálculo de 𝑤𝑥𝑁𝑡 Este cálculo é necessário pois ao se calcular o termo 𝑤𝑥𝑁𝑡1 o método das diferenças finitas centrado requer o conhecimento de 𝑤𝑥𝑁𝑡 e 𝑤𝑥𝑁𝑡 Para isso propõese que o último passo da integração numérica no tempo seja realizado através da aproximação de primeira ordem de Newton Assim temse 𝑤𝑥𝑁𝑡 𝑤𝑥𝑁𝑡1 ℎ𝑡 𝑑𝑤𝑥𝑁𝑡1 𝑑𝑥 39 Para aproximar a primeira derivada da função 𝑤 em função do tempo utilizase a equação 7 o que resulta em 𝑤𝑥𝑁𝑡 𝑤𝑥𝑁𝑡1 ℎ𝑡 𝑤𝑥𝑁𝑡𝑤𝑥𝑁𝑡2 2ℎ𝑡 40 Reorganizando se obtém a relação necessária para a conclusão da integração numérica no tempo 𝑤𝑥𝑁𝑡 2𝑤𝑥𝑁𝑡1 𝑤𝑥𝑁𝑡2 0 41 Devese então organizar os dados e escrever o problema na forma 𝐴𝑤 𝑞 para se obter a solução desejada 14 Caso Dinâmico Viga Em Balanço O procedimento para a solução do caso dinâmico da viga em balanço é parecido com o do caso dinâmico da viga biapoiada O primeiro passo é escrever as condições de contorno com a dependência temporal Isso pode ser feito aproveitando o trabalho que já foi realizado no cálculo das condições de contorno do caso estático resultando nas seguintes relações 𝑤0𝑡 0 𝑤1𝑡 𝑤1𝑡 0 𝑤𝑁𝑥1𝑡 2𝑤𝑁𝑥𝑡 𝑤𝑁𝑥1𝑡 0 42 𝑤𝑁𝑥2𝑡 2𝑤𝑁𝑥1𝑡 2𝑤𝑁𝑥1𝑡 𝑤𝑁𝑥2𝑡 0 Já as condições iniciais podem ser definidas da mesma forma que no problema dinâmico da viga biapoiada Assim escrevese as condições iniciais 𝑤𝑥0 𝑤0𝑥 𝑤𝑥1 𝑤𝑥1 2ℎ𝑡𝑤 0𝑥 43 Já a regra geral do sistema é exatamente a mesma definida na equação 38 Da mesma forma que no caso anterior a conclusão da integração numérica no tempo depende da realização do último passo que pode ser concluída utilizando a relação mostrada na equação 41 17 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Assim no espaço as condições de contorno subsidiam o cálculo dos índices 1 0 𝑁𝑥 1 e 𝑁𝑥 2 Já no tempo as condições iniciais subsidiam o cálculo dos índices 0 e 1 A regra geral da equação 38 trabalha com os índices de 𝑖 2 a 𝑖 2 no espaço e de 𝑗 1 a 𝑗 1 no tempo Assim no espaço a regra geral deve ser aplicada nos índices de 1 a 𝑁𝑥 e no tempo de 1 a 𝑁𝑡 1 Desta forma o número de variáveis a se resolver é de 𝑁𝑥 4𝑁𝑡 2 Novamente o próximo passo é organizar o problema na forma 𝐴𝑤 𝑞 e obter sua solução Em ambos os casos a discretização no espaço e no tempo pode levar a matrizes muito grandes Por isso privilegiase métodos iterativos de s O PROBLEMA DE AUTOVETORES E AUTOVALORES MA DE AUTOVALORES E AUTOVETORES Nosso propósito será determinar certas configurações denominadas autovetores ou modos Buscaremos também as condições críticas para que esses estados sejam possíveis Um parâmetro que descreve uma condição desse tipo é chamado um autovalor Como exemplos existem a determinação de modos de vibração e suas respectivas frequências naturais em sistemas oscilantes e de modos de cargas de flambagem em problemas de estabilidade elástica Problemas Padrão e Generalizado O problema pode ser colocado na forma 𝐴𝑥 𝜆𝑥 Ou seja neste sistema de equações queremos determinar os vetores 𝑥 que são levados pela matriz 𝐴 ao subespaço gerado pelo próprio vetor 𝑥 bem como os valores 𝜆 que satisfazem esta propriedade Para a vasta maioria dos problemas de interesse das engenharias e das ciências exatas a matriz 𝐴 é simétrica e invertível Segue da equação 𝐴𝑥 𝜆𝑥 que 𝐴𝑥 𝜆𝑥 𝐴𝑥 𝜆𝐼𝑥 𝐴𝑥 𝜆𝐼𝑥 0 𝐴 𝜆𝐼𝑥 0 Ou seja tratase de um sistema de equações homogêneo que só possui soluções não triviais se 𝑑𝑒𝑡𝐴 𝜆𝐼 0 A última equação fornece uma equação polinomial de grau 𝑛 em 𝜆 cujas 𝑛 raízes 𝜆𝑖 são os autovalores de𝐴 Se cada um desses valores for substituído na equação 𝐴 𝜆𝐼𝑥 0 então obtémse os 𝑛 vetores 𝑥𝑖 os respectivos autovetores associados aos autovalores 𝜆𝑖 Exemplo Determine os autovalores e os autovetores da matriz 𝐴 responsável pela reflexão de vetores em torno do eixo 𝑥 𝐴 1 0 0 1 Devemos essencialmente trabalhar com a equação dos autovetores 𝐴 𝜆𝐼𝑥 0 e com a equação dos autovalores 𝑑𝑒𝑡𝐴 𝜆𝐼 0 Começamos resolvendo a segunda equação pois precisaremos dos autovalores para substituir na primeira equação a fim de determinarmos os autovetores Vejamos 𝑑𝑒𝑡𝐴 𝜆𝐼 0 Se 𝜆 1 voltamos na equação dos autovetores e temos que Resolvendo este sistema linear temos que a solução é dada por Ou seja encontramos o autovetor 𝑥 1 0 associado ao autovalor 𝜆 1 Com efeito os vetores que são levados para o mesmo subespaço pela matriz 𝐴 são os vetores do eixo x que são gerados logicamente pelo autovetor 𝑥 1 0 18 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Similarmente Se 𝜆 1 voltamos na equação dos autovetores e temos que Resolvendo este sistema linear temos que a solução é dada por Ou seja encontramos o autovetor 𝑥 0 1 associado ao autovalor 𝜆 1 Com efeito os vetores que são levados para o mesmo subespaço pela matriz 𝐴 são os vetores do eixo y que são gerados logicamente pelo autovetor 𝑥 0 1 A equação 𝐴𝑥 𝜆𝑥 é chamada de problema padrão de autovalores e autovetores Uma situação mais completa é o chamado problema generalizado de autovalores e autovetores na forma 𝐴𝑥 𝜆𝐵𝑥 Na qual se admite que as matrizes 𝐴 e 𝐵 são ambas simétricas e pelo menos uma delas é não singular geralmente 𝐵 Considere o seguinte problema 𝐴𝑥 𝜆𝐵𝑥 onde 𝐴 600 600 0 600 1800 1200 0 1200 3000 𝐵 1 0 0 0 15 0 0 0 2 Como 𝐵 é não singular ou seja det𝐵 0 pode se fazer 𝐴𝑥 𝜆𝐵𝑥 𝐵1 𝐴𝑥 𝐵1 𝜆𝐵𝑥 𝐵1𝐴𝑥 𝜆𝐵1𝐵𝑥 𝐵1𝐴𝑥 𝜆𝐼𝑥 𝐵1𝐴𝑥 𝜆𝑥 Ou seja 𝐻𝑥 𝜆𝑥 onde 𝐻 𝐵1𝐴𝑥 Vamos encontrar o autovalor 𝜆 de𝐻 usando o comando spec do scilab Para se obter os autovetores e os autovalores de H fazse A 600 600 0600 1800 12000 1200 3000 B1 0 00 15 00 0 2 HinvBA R diagevalsspecH Saída do Algoritmo Observação Rdiagevals specH retorna na matriz diagonal evals os autovalores e em R os autovetores direitos 19 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Então 𝜆 21087 é o menor autovalor do problema associado ao autovetor 𝑥 08133 05274 02455𝑇 Enquanto 𝜆 212516 é o maior autovalor do problema associado ao autovetor 𝑥 027304 069406 066126𝑇 APLICAÇÕES NAS ENGENHARIAS E CIÊNCIAS A figura abaixo exibe um modelo de um sistema vibratório linear não amortecido de três graus de liberdade com os deslocamentos das massas 𝑢 𝑢1 𝑢2 𝑢3𝑇 Foram adotados os parâmetros 𝑘1 𝑘4 10 𝑁 𝑚 𝑘2 𝑘3 30 𝑁 𝑚 𝑚1 𝑚2 𝑚3 1𝑘𝑔 A energia cinética do sistema é dada por 𝑇 1 2 𝑚1𝑢1 2 𝑚2𝑢2 2 𝑚3𝑢3 2 A energia de deformação é dada por 𝑈 1 2 𝑘1𝑢1 2 𝑘2𝑢2 𝑢12 𝑘3𝑢3 𝑢22 𝑘4𝑢3 2 Aplicando a equação de Lagrange 𝑑 𝑑𝑡 𝑇 𝑢𝑖 𝑈 𝑢𝑖 0 𝑖 123 Chegase a um sistema de três equações diferenciais ordinárias lineares de segunda ordem 𝑴𝒖 𝑲𝒖 𝟎 onde se definem as matrizes de massa e de rigidez Sabemos da teoria de equações diferenciais ordinárias de segunda ordem que a solução de 𝑴𝒖 𝑲𝒖 𝟎 é uma solução periódica da forma 𝒖 𝒖𝒄𝒐𝒔𝒘𝒕 em que 𝒖 é um modo de vibração isto é um vetor cujas componentes variam no tempo mantendo as mesmas relações entre suas componentes Em outras palavras as componentes do modo evoluem no tempo em fase e na mesma frequência 𝒘 em rads Substituindo 𝒖 𝒖𝒄𝒐𝒔𝒘𝒕 na equação diferencial que governa o movimento e cancelando o cosseno temse 𝑲𝒖 𝑤2𝑴𝒖 Tratase um problema de autovalores as frequências ao quadrado e autovetores modos de vibração Colocando na forma padrão temse que 𝑬𝒖 𝑤2𝒖 onde 𝑬 𝑴𝟏𝑲 Para se obter os autovetores e autovalores de 𝑬 fazse Meye33 K 40 30 030 60 300 30 40 EinvMK R diagevals specE O primeiro modo do sistema é o autovetor 𝑥 05543 06207 05543𝑇 associado ao autovalor 𝜆 𝑤2 64110 e portanto associado à sua primeira frequência 𝑤1 253 𝑟𝑎𝑑 𝑠 O terceiro modo do sistema corresponde ao autovetor 04389 0 78403 043891𝑇 associado ao autovalor 𝜆 𝑤2 9358 e portanto associado à sua terceira frequência 𝑤1 9672 𝑟𝑎𝑑 𝑠 Exercícios propostos e aplicações 1 Encontre os polinômios de Taylor 𝑇𝑛 centrados em 𝑎 para 𝑛 2345 utilizando o geogebra para os cálculos das derivadas Então trace esses polinômios e 𝑓 na mesma tela usando o scilab 20 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 a 𝑓𝑥 cot 𝑥 𝑎 𝜋 4 b 𝑓𝑥 1 x2 3 𝑎 0 2 Use a série de Taylor de 𝐿𝑛1 𝑥 para calcular 𝐿𝑛14 com precisão de 5 casas decimais 3 Sob diversas hipóteses simplificadoras a altura estacionária de um lençol freático em um aquífero de águas subterrâneas não confinado unidimensional pode ser modelada pela seguinte EDO 𝐾ℎ 𝑑2ℎ 𝑑𝑥2 𝑁 0 Tal que 𝑥 distância em metros 𝐾 condutividade hidráulica 𝑚𝑑 ℎ altura do lençol freático em metros e ℎ altura média do lençol freático em metros e 𝑁 taxa de infiltração 𝑚𝑑 Determine a altura do lençol freático para 𝑥 0 a 1000 m em que ℎ0 10 m e h10005 m Use os seguintes parâmetros para os seus cálculos 𝐾 1 𝑁 00001 Considere a altura média do lençol freático como sendo a média das condições de contorno Obtenha sua solução com o método numérico das diferenças finitas considerando 𝑥 10 𝑚 4 Resolva numericamente o problema de valor de contorno linear abaixo considerando 𝑛 1 21 pontos e espaçamento h005 usando o método de diferenças finitas centradas 𝑦 2𝑦 𝑦 𝑥 𝑦0 0 𝑦1 1 21 Introdução à Simulação Numérica Lista de Exercícios Teoria e Aplicações 2022 Gabarito dos exercícios propostos 1 2 3 4 x y 0 0 005 01427425 01 02714316 015 03869802 02 04902457 025 05820329 03 06630974 035 07341484 04 07958514 045 08488306 05 08936717 055 0930924 06 09611024 065 09846898 07 10021386 075 1013873 08 102029 085 1021762 09 10186371 095 10112417 1 1