Resolução de uma Edo de 2 a Ordem do Tipo Pvi Via o Método de Euler

Teste Computacional 2 de Algoritmos Numericos Resolucao de uma EDO de 2a ordem do tipo PVI via o metodo de Euler 1 Objetivo O objetivo deste trabalho e usar o metodo de Euler para resolver uma EDO 2a ordem e em particular obter a solucao do problema da deflexao de uma viga engastada 2 O metodo de Euler para um PVI de 2a ordem Um problema matematico bastante comum que aparece aos se modelar fenˆomenos fısicos sao as EDOs de 2a ordem do tipo PVI Problema de Valor Inicial ou seja problemas onde aparece uma equacao diferencial ordinaria de 2a ordem com valores fornecidos em um unico ponto x0 a tal como descrito abaixo y fx y y ya ya ya qa Este PVI pode ser resolvido empregando as mesmas tecnicas que aquelas utilizadas para resolver um sistema de duas equacoes diferencias de 1a ordem Para isto basta reescrever o PVI de 2a ordem como um sistema de equacoes diferenciais ordinarias de 1a ordem fazendo a seguinte transformacao y1 y y2 y pois assim y 1 y y2 y 2 y fx y y Portanto o problema passar a ser obter as funcoes y1x e y2x solucao de um sistema de duas equacoes diferenciais de 1a ordem com valores iniciais definidos abaixo y 1 y2 y 2 fx y1 y2 y1a ya y2a qa Esse sistema pode tambem ser escrito por y 1 f1x y1 y2 y 2 f2x y1 y2 y1a ya y2a qa onde f1x y1 y2 y2 e f2x y1 y2 fx y1 y2 E lembrando que apesar de se obter duas funcoes y1x e y2x o interesse final e apenas saber yx y1x isto e calculase tambem como efeito colaterala funcao y2x que representa a derivada de yx Para resolver um sistema de duas equacoes deve ser usado o metodo de Euler que esta descrito a seguir 21 Sistemas de EDOs de 1a ordem Seja o sistema de duas equacoes diferencias ordinarias de 1a ordem com valores iniciais tal como o definido abaixo y 1 f1x y1 y2 y 2 f2x y1 y2 y1a y1a y2a y2a A solucao analıtica em um domınio D a b e o par de funcoes y1x e y2x que satisfaz as equacoes diferenciais e as condicoes iniciais y1a y1a e y2a y2a em D a b Existem varios metodos para obter a solucao de forma numerica Metodos conhecidos como metodos de Runge Kutta sao bastante empregados na pratica Neste trabalho sera usado apenas um dos metodos mais simples o metodo de Euler que e tambem um metodo Runge Kutta de 1a ordem O metodo de Euler para duas equacoes Quando ha duas equacoes as aproximacoes numericas para y1x e y2x sao calculadas sincronizadamente para y1x e y2x Nos metodos de passo simples a cada passo da discretizacao do domınio calcula se as aproximacoes para y1xi 1 e para y2xi 1 ou seja os valores y1i1 e y2i1 a partir dos valores numericos em xi isto e a partir de y1i e y2i O metodo de Euler e dado por y1i1 y1i hf1xi y1i y2i y2i1 y2i hf2xi y1i y2i Que pode tambem ser escrito como y1i 1 y1i h f1 y2i 1 y2i h f2 onde f1 f1xi y1i y2i e a inclinacao da funcao y1 em xi f2 f2xi y1i y2i e a inclinacao da funcao y2 em xi h b am e m e a particao do dominio 3 Implementacoes execucoes e problemas a serem tratados 5 pts Usando como modelo uma possıvel implementacao em python do metodo Euler resolver para sistemasde EDOs fornecida ou seja usando o codigo EulerSistemas2Eqpy resolver os problemas descritos abaixo 1 Resolucao de um problema de validacao problema A Determinar a solucao de y 1 f1x y1 y2 y2 y 2 f2x y1 y2 2y1 y2 x2 y100 10 y200 00 em D 0 1 com m 10 subintervalos Se vocˆe quiser fazer alguma verificacao as solucoes esperadas exatas sao y1exata 025e2x 2x2 2x 3 y2exata y 1 2 Calculo da deflexao de uma viga engastada problema B Seja uma viga de comprimento L engastada em uma das suas extremidades Seja uma carga pontual P sendo aplicada em um ponto que esta a x unidades de medida do ponto onde a viga esta presa A deflexao na extremidade livre x L depende de onde a carga esta sendo aplicada Quando a carga e aplicada bem proximo a extremidade fixa a deflexao e menor do que quando e aplicada mais proxima a extremidade livre e esta deflexao aumenta a medida que a carga se distancia da borda engastada Vale observar que quando a carga e aplicada no ponto onde a viga esta presa em x 0 a deflexao e nula ou seja y0 0 Denominando de y a deflexao na extremidade livre x L o modelo matematico que descreve esta deflexao e dado por y P EI 1 y215L x y0 0 y0 0 onde E e o modulo de Young do material e I o momento de inercia da secao da viga Obter a deflexao de uma viga engastada em aco de tamanho L 120cm sujeita a acao de uma carga pontual P 15N aplicada em um ponto x Considere que o modulo de Young do aco e E 200000Ncm2 e que o momento de inercia da secao da viga vale 25cm2 Observe que o PVI de 2a ordem que descreve o fenˆomeno y P EI 1 y215L x y0 0 y0 0 pode ser reescrito por um sistema de duas equacoes diferenciais de 1a ordem fazendo a seguinte transformacao y1 y y2 y com isso o PVI passa a ser um sistema de duas equacoes dado abaixo y 1 y y2 y 2 y P EI 1 y2215L x y1x 0 0 y2x 0 0 Obter a solucao usando uma particao com m 40 Obter o valor da deflexao o valor de y1 y para todos os pontos da particao ou seja o valor de y1 y em funcao de x 4 Menu de entrada e valores de saıda Seu programa deve exibir na tela um menu para o usuario similar ao mostrado abaixo Digite uma opcao A Resolver o problema A exemplo da validacao B Resolver o problema B determincao da deflecao na viga S Sair Escolha Apos a escolha da opcao o usuario deve fornecer os dados necessarios para a resolucao do problema ou seja a b y1a y2a e m Fornecidos os dados de entrada o codigo deve exibir na tela a solucao numerica do problema escolhido com a particao estabelecida incluindo o ponto x a isto e para cada xi deve se exibir os valores de y1i e y2i Por exemplo Suponha um domınio D 10 20 com m 2 subdivisoes e valores fictıcios para y1 e y2 A saıda seria x 10 y1 40 y2 07 x 15 y1 385 y2 06 x 20 y1 373 y2 05 O formato de saıda fica a seu criterio mas o codigo deve ter como saıda os valores de y1 e y2 para todos os pontos da particao 5 Relatorio 5 pts Alem de me enviar o codigo em python devese relatar os resultados obtidos apresentando os problemas que estao sendo resolvidos Para fazer este relatorio siga aproximadamente os seguintes passos a Apresentar a solucao numerica do problema A Exibir os valores y1i e y2i obtidos com m 10 subintervalos b Apresentar o problema B isto e o problema da deflexao da viga e exibir os valores de yi para m 40 6 Condicoes de entrega 61 Grupo Este trabalho devera ser realizado em grupos de no maximo 3 alunos OBS Escreva os nomes dos componentes no codigo principal do seu programa Introdução A resolução de sistemas de equações diferenciais ordinárias EDOs é uma parte essencial da modelagem em diversas áreas como física biologia e engenharia Dentre os métodos numé ricos disponíveis para resolver esses sistemas o Método de Euler se destaca pela sua simplici dade e eficiência Ele é frequentemente empregado para aproximar soluções de EDOs quando a solução exata é desconhecida ou difícil de calcular Para sistemas de EDOs de primeira ordem esse método permite calcular as aproximações numéricas das variáveis dependentes de forma simultânea A fórmula do método de Euler para um sistema de duas equações diferenciais ordinárias por exemplo pode ser expressa da seguinte maneira y1i1 y1i hf1xiy1iy2i y2i1 y2i hf2xiy1iy2i onde f1xiy1iy2i é a inclinação da função y1 em xi f2xiy1iy2i é a inclinação da função y2 em xi e h é o tamanho do passo de integração Neste trabalho dois sistemas de equações diferenciais de ordem 2 serão resolvidos utili zando o método de Euler As resoluções serão obtidas a partir de código implementado em linguagem Python que aplica o método numérico escolhido para resolver ambos os problemas Conforme proposta da atividade as seções a seguir tratam da resolução numérica de dois pro blemas dados na forma de sistemas de equações diferenciais de primeira ordem O primeiro problema apresenta uma vertente experimental consistindo na ideia de verificar se o código implementado e o próprio método em si apresenta boa funcionalidade o segundo problema por sua vez trata de um exemplo prático com parâmetros físicos definidos na descrição das instruções da atividade 1 Resolução de um problema de validação problema A Neste problema desejase resolver o seguinte sistema de equações diferenciais ordinárias y 1 f1xy1y2 y2 y 2 f2xy1y2 2y1 y2 x2 y100 10 y200 00 no intervalo 01 com m 10 subintervalos ou partições O sistema será resolvido numericamente por meio do Método de Euler assistido por có digo construído em Python com configuração disposta no Apêndice Para avaliar a acurácia do método a solução numérica será comparada com a solução exata dada por y1x 025e2x 2x2 2x3 y2x 05e2x x05 Tal comparação consistiu simplesmente do cálculo do erro absoluto local entre soluções numérica e exata e permite verificar a precisão da solução obtida numericamente Para atender às especificações desta parte do trabalho definimos a solução exata do trabalho e as funções f1 e f2 que descreve o sistema O script 1 apresenta essas definições voltadas exclusivamente para o problema A 1 Problema A 2 Solucao exata do problema A 3 def y1exatax 4 return 025mathexp2x 2x2 2x 3 5 def y2exatax 6 return 05mathexp2x x 05 7 8 Funcoes do problema A 9 def f1Ax y1 y2 10 return y2 11 def f2Ax y1 y2 12 return 2y1 y2 x2 Script 1 Parte do código que define a solução exata para o problema sob as condições especi ficadas e as funções que compõem o sistema A etapa seguinte consite na impressão de um menu inicial que apresenta as opções de resolução do problema A ou do problema B será abordado na próxima seção ou sair da execução do programa Para impressão deste menu utilizouse a parte do código expressa no Script 2 1 print Escolha uma opcao 2 printA Resolver o problema A exemplo de validacao 3 printB Resolver o problema B determinacao da deflecao na viga 4 printS Sair 5 opcao input stripupper 6 print Script 2 Apresentação do menu de opções disponíveis para execução Essas opções são impressas dentre de um loop while para garantir que a saída ocorra ape nas quando o usuário escolher a opção Sair Além disso utilizamos o comando stripupper para remove os caracteres em branco do inicio e do final da string e transforma todos os carac teres em maiusculas para tornar válida a escolha das letras correspondentes na forma minúscula e aumentar assim o número de escolhas válidas No caso de a opção A ser escolhida solicitase o fornecimento dos dados necessários para a solução do problema a b y1a y2a e m De posse dessas informações o código aplica o método de Euler e apresenta as estimativas obtidas impressas em tela Para fins de exemplificação considerando as entradas a 0 b 1 y1a 1 y2a 0 e m 10 obtemos as informações dispostas na Tabela 1 Tabela 1 Resultados obtidos para o problema A quando a 0 b 1 y1a 1 y2a 0 e m 10 xi y1i y2i y1xiy1i y2xiy2i 000 10000 00000 0 0 010 10000 02000 00104 00107 020 10200 04190 00230 00269 030 10619 06609 00386 00502 040 11280 09304 00584 00824 050 12210 12330 00835 01261 060 13443 15755 01157 01845 070 15019 19659 01569 02617 080 16985 24139 02098 03626 090 19399 29310 02776 04938 100 22330 35310 03643 06635 Como a primeira linha apresenta os valores iniciais do problema é óbvio que não há dife rença entre o resultado numérico e exato nessa posição Nas demais posições são percebidos erros que vão aumentando à medida em que avançamos em posições mais distântes da posi ção inicial Este comportamento é esperado pois os erros produzidos vão sendo acumulados no processo de resolução uma vez que partimos de uma estimativa que apresenta erros para estimar um novo valor em posição distinta Convém destacar que a solução exata incluída no código é válida somente para as condições iniciais apresentadas na Tabela 1 Caso as condições iniciais sejam diferentes a solução exata considerada não tem validade comparativa com essa configuração 2 Cálculo da deflexão de uma viga engastada problema B O problema B trata da deflexão de uma viga engastada sujeita a uma carga pontual Este problema é uma aplicação direta de equações diferenciais na engenharia estrutural especial mente no cálculo de deformações em vigas A viga considerada no problema possui compri mento L 120 cm e está sujeita a uma carga pontual P 15 N As características do material da viga que é de aço incluem o módulo de Young E 200000 Ncm² e o momento de inércia da seção transversal I 25 cm² O sistema de equações diferenciais ordinárias de ordem que descreve a deflexão y1 e a inclinação y2 da viga obtidos pela simplicação do modelo primário que é diferencial de ordem superior é dado por y 1 f1xy1y2 y2 y 2 f2xy1y2 P EI1y2215Lx y100 00 y200 00 A equação diferencial acima representa a forma geométrica da viga sob a ação da carga O termo 1 y2 215 é uma não linearidade que surge devido à inclinação da viga e que pode tornar a resolução analítica do problema uma tarefa complexa Dentre do código os valores dos parâmetros físicos estão definidos conforme Script 3 mas podem ser modificados para outros parâmetros a depender do desejo ou necessidade de quem usa o código 1 Definindo os parametros do problema B viga de aco 2 L 120 cm 3 P 15 N 4 E 200000 Ncm modulo de Young do aco 5 I 25 cm4 momento de inercia Script 3 Parâmetros físicos e estruturais considerados para o problema B Assim como no problema anterior para resolver este sistema numericamente utilizamos o método de Euler O método discretiza o domínio e aproxima as soluções de y1 deflexão e y2 inclinação em cada ponto xi A resolução do Problema B ilustra a aplicação prática de EDOs em problemas de enge nharia especialmente na análise de estruturas sujeitas a cargas O Método de Euler apesar de sua simplicidade oferece estimativas que podem ser úteis em contextos práticos e como dados iniciais para métodos mais avançados e mais precisos Para fins de exemplificação de funcionamento do código considerouse as entradas a 0 b 120 y1a 0 y2a 0 e m 40 obtemos as informações dispostas na Tabela 1 Os resultados apresentados de fato corroboram com os comportamentos físicos esperados pois o valor da deflexão aumenta à medida em que avançamos no comprimento da viga Observe também que o problema tem potencial para ser expandido para componentes diferenciados e diversificados Inclusive o código pode ser facilmente modificado para atender problemas de Tabela 2 Deflexões obtidas para o problema B quando a 0 b 120 y1a 0 y2a 0 e m 40 xi y1i xi y1i 000 00000 6300 05727 300 00000 6600 06237 600 00032 6900 06763 900 00096 7200 07304 1200 00191 7500 07858 1500 00316 7800 08425 1800 00470 8100 09004 2100 00652 8400 09595 2400 00862 8700 10196 2700 01098 9000 10807 3000 01361 9300 11426 3300 01648 9600 12054 3600 01960 9900 12689 3900 02296 10200 13331 4200 02654 10500 13978 4500 03034 10800 14631 4800 03435 11100 15287 5100 03856 11400 15946 5400 04296 11700 16608 5700 04756 12000 17272 6000 05233 ordens maiores apenas com a inserção de funções que representa outras equações diferenciais de primeira ordem obtidas por redução de ordens e o incremento de poucas linhas na função que aplica o método de Euler Uma das vantagens do uso desses métodos extraídas da experiência com esta resolução de problemas é a facilidade de estimar soluções para problemas distintos alterando poucos para mêtros nos códigos funções que definem o problema condições iniciais tamanho do passo Concluímos que soluções numéricas para equações diferenciais são muito úteis e importantes para resolver problemas em diversas áreas do conhecimento Apêndice 1 import math 2 3 Mensagem inicial 4 def mensageminicial 5 print Iniciando o codigo 6 7 Problema A 8 Solucao exata do problema A 9 def y1exatax 10 return 025mathexp2x 2x2 2x 3 11 def y2exatax 12 return 05mathexp2x x 05 13 14 Funcoes do problema A 15 def f1Ax y1 y2 16 return y2 17 def f2Ax y1 y2 18 return 2y1 y2 x2 19 20 Problema B 21 Funcoes do problema B deflexao de uma viga 22 def f1Bx y1 y2 23 return y2 24 def f2Bx y1 y2 P EI L 25 return PEI1 y2215L x 26 27 Esquema Numerico 28 Funcao que aplica o metodo de Euler para sistemas de EDOs 29 def eulera b y1a y2a m f1 f2 paramsNone compararFalse 30 printObtendo a solucao via o Metodo de Euler 31 y1 y1a Condicao inicial para y1 32 y2 y2a Condicao inicial para y2 33 h b am Tamanho do passo 34 x a Posicao de partida 35 36 printf Dados Para x x2f y1 y1a4f y2 y2a4 f condicao inicial 37 38 Aplica o esquema numerico do metodo de Euler 39 for i in range1 m 1 40 if params 41 y1viz y1 hf1x y1 y2 42 y2viz y2 hf2x y1 y2 params 43 else 44 y1viz y1 hf1x y1 y2 45 y2viz y2 hf2x y1 y2 46 47 Atualizacao dos valores 48 x x h 49 y1 y1viz 50 y2 y2viz 51 52 if comparar 53 y1ex y1exatax 54 y2ex y2exatax 55 printfPasso i Para x x2f y1 y14f y2 y24f dy1 absy1 y1ex4f dy2 absy2 y2ex4f 56 else 57 printfPasso i Para x x2f y1 y14f y2 y24f 58 59 return y1 y2 60 61 Funcao principal 62 def menu 63 mensageminicial 64 65 while True 66 print Escolha uma opcao 67 printA Resolver o problema A exemplo de validacao 68 printB Resolver o problema B determinacao da deflecao na viga 69 printS Sair 70 opcao input stripupper 71 stripupper remove os caracteres em branco do inicio e do final da string e transforma todos os caracteres em maiusculas para facilidade de validacao da escolha 72 print 73 74 if opcaoA 75 Entrada do usuario para o problema A 76 a floatinputDigite o limite inferior do intervalo a 77 b floatinputDigite o limite superior do intervalo b 78 y1a floatinputDigite o valor de y1a 79 y2a floatinputDigite o valor de y2a 80 m intinputDigite a quantidade de subdivisoes m 81 82 printf Resolvendo problema A no intervalo a b com m m subintervalos 83 if y1a1 and y2a0 84 eulera b y1a y2a m f1A f2A compararTrue 85 else 86 eulera b y1a y2a m f1A f2A compararFalse 87 88 elif opcao B 89 Definindo os parametros do problema B viga de aco 90 L 120 cm 91 P 15 N 92 E 200000 Ncm2 modulo de Young do aco 93 I 25 cm4 momento de inercia 94 EI EI Produto EI 95 96 Entrada do usuario para o problema B 97 a floatinputDigite o limite inferior do intervalo a 98 b floatinputDigite o limite superior do intervalo b 99 y1a floatinputDigite o valor de y1a 100 y2a floatinputDigite o valor de y2a 101 m intinputDigite a quantidade de subdivisoes m 102 103 printfResolvendo problema B no intervalo a b com m m subintervalos 104 eulera b y1a y2a m f1B f2B paramsP EI L compararFalse 105 106 elif opcao S 107 printSaindo 108 break 109 110 else 111 printOpcao invalida 112 113 Inicia o menu 114 menu 115 printFIM