Lista - Algoritmos Numéricos 2023 2

|f(x_3)| = | - 3.331278 | \leq 0.001 (Falso!). A condição é falsa, por isso o processo continua \left(x_4 = x_3-\frac{(x_3-x_2)f(x_3)}{f(x_3)-f(x_2)}\right) até ambas as condições do critério de paragem se verificarem simultaneamente, que acontece no final da 3ª iteração. Assim, x_4 = 0.117575, f(x_4) = 0.059767, x_5 = 0.117609 e f(x_5) = -0.000255. A verificação do critério de paragem será: \left|\frac{x_5-x_4}{x_5}\right| = \frac{|0.117609-0.117575|}{|0.117609|} = 0.000292 < 0.001 (Verdadeiro!) |f(x_5)| = 0.000255 < 0.001 (Verdadeiro!) A solução é encontrada ao fim de 3 iterações: A^* \approx 0.117609 radianos. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = pi*pot((300/cos(x1)),2)*0.8/(0.5*pi*14*14*(1+sen(x1)-0.5*cos(x1)))-1200 Valores Iniciais x1 = 0.000000 x2 = 0.125664 Tolerâncias do Critério de Paragem e1 = 1.000000e-03 e2 = 1.000000e-03 Número Máximo de Iterações NMAX = 4 Iteração 1 xn = 0.119521 fxn = -3.331278 ... Iteração 3 xn = 0.117609 fxn = -0.000252 Número de Iterações Realizadas = 3 Solução xn = 0.117609 fxn = -0.000252 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR 2. Um certo equipamento de 20000 euros vai ser pago durante 6 anos. O pagamento anual é de 4000 euros. A relação entre o custo do equipamento P, o pagamento anual A, o número de anos n e a taxa de juro i é a seguinte: A = \dfrac{P\ \ i(1+i)^n}{\ (1+i)^n-1} . Utilize o método iterativo mais adequado para determinar a taxa de juro utilizada nos cálculos. O valor da taxa de juro pertence ao intervalo [0.05, 0.15]. Para a paragem do processo iterativo use \varepsilon_1 = \varepsilon_2 = 0.05 ou no máximo 3 iterações. Resolução: Mudança de variável: i \rightarrow x e A \rightarrow f. Utiliza-se o Método da Secante, por simplificar os cálculos (não utiliza derivadas) e por isso usam-se dois valores iniciais, x1 = 0.05 e x2 = 0.15. Coloca-se a função na forma f(x) = 0, i.e., f(x) = 20000\ \dfrac{x(1+x)^6}{(1+x)^6-1} - 4000 = 0 Realiza-se a primeira iteração do método, utilizando (A.3), com k = 2: x_3 = x_2 - \dfrac{(x_2-x_1)f(x_2)}{f(x_2)-f(x_1)} = 0.15 - \dfrac{(0.15-0.05)\times 1284.738131}{1284.738131 - (-59.650638)} = 0.054437 Usa-se (A.4) para verificar o critério de paragem. |f(x_3)| = | - 3.562855 | \leq 0.05 (Falso!) Não se verificou a segunda condição do critério de paragem e por isso não se calcula a primeira \left|\dfrac{x_3-x_2}{x_3}\right| \leq \varepsilon_1\right), pelo que as iterações continuarão a realizar-se até que ambas as condições sejam satisfeitas, que acontece no final da 3ª iteração. 4) (valor = 3 pontos) Aplique o método das diferenças finitas para resolver numericamente o seguinte problema de equações diferenciais com condições de contorno, aplicando uma grade de 7x7 pontos internos com os respectivos valores (x,y): (1.000, 0.000) (1.000, 0.125) (1.000, 0.250) (1.000, 0.375) (1.000, 0.500) (1.000, 0.625) (1.000, 0.750) (1.000, 0.875) (1.000, 1.000) (0.875, 0.000) (0.875, 0.125) (0.875, 0.250) (0.750, 0.000) (0.750, 0.125) (0.750, 0.250) (0.625, 0.000) (0.625, 0.125) (0.625, 0.250) (0.500, 0.000) (0.500, 0.125) (0.500, 0.250) (0.375, 0.000) (0.375, 0.125) (0.375, 0.250) (0.250, 0.000) (0.250, 0.125) (0.250, 0.250) (0.125, 0.000) (0.125, 0.125) (0.125, 0.250) (0.000, 0.000) (0.000, 0.125) (0.000, 0.250) (0.000, 0.000) (0.000, 0.125) (0.000, 0.250) (0.000, 0.375) (0.000, 0.500) (0.000, 0.625) Será gerado um sistema de equações lineares. A solução deste poderá ser feita através do python ou do excel. OBS: se você apenas indicar claramente como serão calculados os valores dos coeficientes e do lado direito do sistema de equações lineares resultante (sem calcular os respectivos valores, efetivamente, e – portanto – sem determinar os valores de U(x,y)), será atribuído 50% do valor da questão. Esta resposta poderá ser dada em termos matemáticos ou através de um programa python ou C++. se você indicar como são calculados os coeficientes e gerar os valores dos mesmos, será atribuído 70% do valor da questão; se você indicar como são gerados os coeficientes, calcular os seus valores e resolver o sistema de equações lineares, será atribuído 100% do valor da questão. A equação diferencial parcial, propriamente dita, é a que se segue: (∂^3 U/∂x^3) + (∂^3 U/∂y^3) + (∂^2 U/∂x^2) + (∂^2 U/∂y^2) = 0 U(0, y) = y^3 if 0 <= y <= 1, U(1, y) = 1 + y^3 - 3y - 3y^2 if 0 <= y <= 1, U(x, 0) = x^3 if 0 <= x <= 1, U(x, 1) = 1 + x^3 - 3x^2 - 3x if 0 <= x <= 1. Solução analítica: U(x, y) = x^3 + y^3 - 3x^2y - 3xy^2. Você pode optar entre usar aproximação e baixa ou de alta precisão para as derivadas numéricas. Não fará diferença na avaliação da nota. 4 Métodos Numéricos: exercícios resolvidos aplicados à Engenharia e outras Ciências Maria Teresa Torres Monteiro com a colaboração de Sara Tribuzi M. N. Morais Universidade do Minho Fevereiro 2012 5) (valor = 3 pontos) Aplique o método das diferenças finitas para resolver numericamente o seguinte problema de equações diferenciais com condições de contorno, aplicando uma grade de 7x7 pontos internos (como no enunciado da questão acima). Será gerada um sistema de equações lineares cuja solução poderá ser feita através do python, octave ou do excel. OBS: • se você apenas indicar claramente como serão calculados os valores dos coeficientes e do lado direito do sistema de equações lineares resultante (sem calcular os respectivos valores, efetivamente, e – portanto – sem determinar os valores de U(x,y)), será atribuído 50% do valor da questão. ▪ Esta resposta poderá ser dada em termos matemáticos ou através de um programa python ou C++. • se você indicar como são calculados os coeficientes e gerar os valores dos mesmos, será atribuído 70% do valor da questão; • se você indicar como são gerados os coeficientes, calcular os seus valores e resolver o sistema de equações lineares, será atribuído 100% do valor da questão. 3 Aos meu alunos, de ontem, de hoje, de sempre... Longo é o caminho ensinado pela teoria. Curto e eficaz, o do exemplo. Séneca. Não se aprende bem a não ser pela própria experiência. Bacon. Prefácio Os Métodos Numéricos têm um papel estrutural, de carácter transversal na formação em cursos de Engenharia. São ainda muito utilizados em problemas nas áreas de Economia, Medicina, Física, Química, entre outras. Os Métodos Numéricos procuram desenvolver processos de cálculo (algoritmos), utilizando uma sequência finita de operações aritméti- cas básicas, de forma a que certos problemas matemáticos se tornem exequíveis. Estes algoritmos envolvem, em geral, um grande número de cálculos aritméticos. Não é pois de estranhar que, nas últimas décadas, com o rápido crescimento das potencialidades dos computadores digitais, o papel dos Métodos Numéricos na resolução de problemas de en- genharia tenha sofrido grande incremento. As mais de duas décadas a leccionar a disciplina de Métodos Numéricos a cursos de Engenharia, levou à constatação da dificuldade dos alunos, em geral, na identificação e utilização/aplicação dos vários tópicos de Métodos Numéricos nas várias unidades curri- culares dos respectivos cursos e ainda noutras situações. Este trabalho surge então, da solicitação ao longo dos anos, pelos alunos, de um documento deste teor, com resolução simples e directa de exercícios de Métodos Numéricos aplicados à Engenharia e outras Ci- ências. É uma publicação de carácter pedagógico que compila uma colecção de exercícios de Métodos Numéricos, simples, resolvidos e todos com aplicação prática. A ideia con- siste na aplicação de métodos numéricos em casos práticos simples, capacitando os alunos para a sua utilização em situações futuras mais complexas. Pela clareza e simplicidade na resolução dos problemas, pode ser utilizado como complemento a um estudo teórico, permitindo uma melhoria no desempenho na unidade curricular de Métodos Numéricos, bem como uma ajuda noutras unidades curriculares e em problemas práticos em geral. i ii PREFÁCIO Características Está estruturada em oito capítulos e um formulário de apoio. São compilados 67 exercí- cios sobre várias temáticas, em que todos eles incluem aplicação prática, apresentando um nível de complexidade compatível com a sua resolução em contexto de sala de aula/exame para um nível de estudos de 1o ciclo. Em cada exercício é incluída uma figura para uma melhor interacção com o utilizador e para conferir uma menor austeridade ao documento. É incluída uma resolução manual de todos os exercícios, bem como uma resolução utilizando o software CONUM e/ou MATLAB. Para simplificação da resolução dos problemas, e uma vez que este documento tem um carácter não-teórico, foi incluído em anexo um formulário de apoio para todo o conteúdo programático abordado no mesmo. São incluídas chamadas oportunas ao formulário ao longo da resolução dos exercícios. Das variadas temáticas existentes em Métodos Numéricos, foram seleccionadas as mais relevantes para a aprendizagem dos alunos de Engenharia ao longo da unidade curricular de Métodos Numéricos. Os oito temas abordados foram divididos em capítulos: • Capítulo 1 - Erro Numérico - Fórmula fundamental do erro • Capítulo 2 - Solução de uma equação não-linear - Métodos iterativos de Newton e da Secante • Capítulo 3 - Sistemas de equações lineares - Método directo da Eliminação de Gauss com Pivotagem Parcial (EGPP) - Método iterativo de Gauss-Seidel • Capítulo 4 - Sistemas de equações não-lineares - Método iterativo de Newton • Capítulo 5 - Interpolação Numérica iii - Polinómio interpolador de Newton (ou das diferenças divididas) - Interpolação segmentada - Splines cúbicas • Capítulo 6 - Mínimos quadrados - Modelos lineares polinomiais e não-polinomiais • Capítulo 7 - Integração numérica - Fórmulas simples e compostas de Simpson, Trapézio e três oitavos • Capítulo 8 - Diferenciação Numérica - Equações diferenciais ordinárias (com condições iniciais e de fronteira) Torna-se completamente absurdo, num contexto de problemas de engenharia, dissociar os Métodos Numéricos do computador. Desta forma, paralelamente à resolução manual dos exercícios propostos, foi utilizado software numérico para devolver resultados compu- tacionais dos problemas (CONUM e MATLAB). CONUM (COmputação NUMérica) é uma aplicação escrita na linguagem C++ para ambiente Windows, desenvolvida pela Professora Edite M. G. P. Fernandes e por José Filipe S. R. Soares com o apoio da Universidade do Minho, no ano lectivo 1993/1994. O software é gratuito e está disponível para download em [1], bem como o livro COmputação NUMérica [3], da mesma autora. MATLAB [2] é uma linguagem de alto-nível e ambiente interactivo que permite realizar tarefas computacionalmente intensivas mais rapidamente que as linguagens de programação tradicionais, como C, C++ ou Fortran. A precisão numérica foi um dos parâmetros com muita influência neste projecto. Os exercícios não foram resolvidos manualmente sempre com a mesma precisão numérica - uns com 6 casas decimais, outros com 4 ou 5. A diferença entre as casas decimais utilizadas na resolução manual e as do software utilizado foi notória na resolução de alguns exercícios. Mesmo entre CONUM e MATLAB, verificaram-se diferenças de precisão. iv PREFÁCIO Por último, a problemática dos Métodos Numéricos utilizados em cada software teve efeito nas resoluções, já que cada um pode utilizar diferentes variantes de um determinado método, não correspondendo por vezes aos utilizados nas resoluções manuais. Agradecimentos Este trabalho enquadra-se num projecto do Grupo de Investigação NSOS (Nonlinear Systems Optimization and Statistics group) do Centro de Investigação Algoritmi da Uni- versidade do Minho para uma Bolsa de Integração na Investigação (BII) financiada pela Fundação para a Ciência e Tecnologia (FCT). O primeiro agradecimento é dirigido à FCT pela oportunidade concedida. Esta oportunidade permitiu à Sara Tribuzi de Morais (bol- seira BII) a integração neste projecto - a sua colaboração foi imprescindível para a concre- tização deste trabalho. Agradece-se também ao Grupo de Métodos Numéricos e Estatísticos do Departamento de Produção e Sistemas, do qual faço parte, alguns dos enunciados dos exercícios apresenta- dos. Ao longo destes anos tivemos sempre a preocupação de inovar, procurando exercícios diferentes, com aplicação prática, para colocar em provas de avaliação e nas fichas de exercícios de trabalho. Nota aos utilizadores Sendo esta a primeira edição impressa, agradece-se que se for detectada alguma incor- recção seja enviada uma mensagem para tm@dps.uminho.pt. Conteúdo Prefácio i 1 Erro numérico 1 2 Solução de uma equação não linear 5 3 Sistemas de equações lineares 35 4 Sistemas de equações não lineares 67 5 Interpolação Numérica 87 6 Mínimos Quadrados 111 7 Integração Numérica 125 8 Diferenciação Numérica 149 A Formulário 185 v Capítulo 1 Erro numérico Neste capítulo pretende-se sensibilizar os alunos para os erros que surgem devido à repre- sentação dos números no computador ou na calculadora e os erros resultantes das operações numéricas efectuadas. Os dados experimentais, obtidos a partir de aparelhos de medição falíveis, vêm afectados de erro. Pretende-se averiguar de que forma é que esse erro vai afectar o resultado de qualquer operação efectuada com esses dados. A fórmula fundamental dos erros (A.1) permite calcular limites superiores dos erros absolutos e relativos das operações aritméticas que utilizam esses dados. Vão ser resolvidos exercícios simples, apenas manualmente, utilizando a fórmula fundamental dos erros. 1 2 CAPÍTULO 1. ERRO NUMÉRICO 1. Calcule um limite superior do erro absoluto e do erro relativo no cálculo da expressão f(x, y, z) = −x + y2 + sen(z), sabendo que são usados os seguintes valores aproximados: x = 1.1 (δx = 0.05); y = 2.04 (δy = 0.005); z = 0.5 rad. (δz = 0.05). Quantos algarismos significativos tem o valor calculado de f? Resolução: A máquina deve ser colocada em modo Radianos. Valor calculado f = 3.5410255386; δx = 0.05, δy = 0.005, δz = 0.05. I =          1.05 ≤ x ≤ 1.15 (intervalo de incerteza para x) 2.035 ≤ y ≤ 2.045 (intervalo de incerteza para y) 0.45 ≤ z ≤ 0.55 (intervalo de incerteza para z) Cálculo dos majorantes: δf δx = −1, |δf δx |I = Mx = 1, δf δy = 2y, |δf δy |I = My = 4.09, δf δz = cos(z), |δf δz |I = Mz = 0.90044712 Fórmula fundamental do erro (A.1): δf ≤ 1 × (0.05) + 4.09 × (0.005) + 0.90044712 × (0.05) = 0.115472355 Limites superiores dos erros absoluto/relativo: δf ≤ 0.115472355, δf |f| ≤ 0.032609. Para identifição dos algarismos significativos deve encontrar-se o primeiro valor su- perior a δf na forma 0.5× potências de base 10. Neste caso tem-se 0.115472355 ≤ 0.5 × 100. Então, colocando o valor de f e o majorante de δf em termos de 0.5× potências de base 10, com o mesmo expoente vem, 3. 5410255386 ×100 0. 5 ×100 Conclui-se que apenas existe 1 algarismo significativo (3) - os algarismos de f cuja posição está à esquerda do algarismo 5 de δf. 3 2. Uma corrente eléctrica atravessa uma resistência R de 20Ω. A resistência foi medida com um erro relativo que não excede 0.01. A intensidade da corrente I é 3.00 ± 0.01 A. Sabendo que a tensão da corrente é dada por V = RI, determine um limite superior do erro absoluto no cálculo da tensão da corrente. Quantos algarismos significativos garante para o valor calculado da tensão? Resolução: Função: V = RI Fórmula fundamental do erro (A.1): δV ≤ MRδR + MIδI Cálculo das derivadas parciais: ∂V ∂R = I, ∂V ∂I = R R = 20Ω; δR R ≤ 0.01 ⇔ δR ≤ 0.2 (limite superior do erro relativo em R). I = 3.00 ± 0.01 ⇒ δI ≤ 0.01 (limite superior do erro absoluto em I). Intervalo de incerteza:    20 − δR ≤ R ≤ 20 + δR ⇔ 19.8 ≤ R ≤ 20.2 (intervalo de incerteza para R) 3 − δI ≤ I ≤ 3 + δI ⇔ 2.99 ≤ I ≤ 3.01 (intervalo de incerteza para I) Cálculo dos majorantes MR, MI no intervalo de incerteza: |∂V ∂R| ≤ MR ⇒ MR = 3.01; |∂V ∂I | ≤ MI ⇒ MI = 20.2 Substituição na fórmula: δV ≤ 20.2 × 0.01 + 3.01 × 0.2 = 0.804 × 100 ≤ 0.5 × 101 V = 20 × 3 = 60 = 6.0 × 101 Conclui-se que apenas existe 1 algarismo significativo - o 6. 4 CAPÍTULO 1. ERRO NUMÉRICO 3. Pretende-se calcular a área de um terreno circular, de raio aproximadamente igual a 250m. Usando 3.14 para valor aproximado de π, quan- tos algarismos significativos apresenta o valor da área? Resolução: Área como função de π e r: f(π, r) = πr2. Derivadas parciais: ∂f ∂π = r2, ∂f ∂r = 2πr δπ ≤ 0.005 δr ≤ 0.5    250 − 0.5 ≤ r ≤ 250 + 0.5 ⇔ 249.5 ≤ r ≤ 250.5 (intervalo de incerteza para r) 3.14 − 0.005 ≤ π ≤ 3.14 + 0.005 ⇔ 3.135 ≤ π ≤ 3.145 (intervalo de incerteza para π) Majorantes: Mπ = 250.52 = 62750.25, Mr = 2 × (3.145) × 250.5 = 1575.645. Fórmula fundamental do erro (A.1): δf ≤ δπMπ + δrMr δf ≤ 0.005×62750.25+0.5×1575.645 = 1101.57375 = 0.110157375×104 ≤ 0.5×104 O valor da área é 196250 = 19.6250 × 104. 19. 6250 ×104 0. 5 ×104 No valor da área (196250) apenas são algarismos significativos o 1 e o 9 (estão na posição à esquerda do 5 de δf). Capítulo 2 Solução de uma equação não linear Neste capítulo vão resolver-se exercícios para obter as soluções de equações algébricas não lineares. O objectivo é determinar os valores de x ∈ ℜ que anulam a função f, i.e., resolver a equação f(x) = 0. Encontrar o conjunto de soluções da equação anterior não é tarefa fácil e não pode ser obtido por métodos algébricos. A solução deste tipo de problema é obtida graças a mé- todos numéricos iterativos. Vão ser utilizados os métodos de Newton (A.2) e da Secante (A.3). Nos processos iterativos a solução exacta só é encontrada após um número infinito de iterações. O processo iterativo é interrompido através da satisfação do critério de pa- ragem (A.4), que deve ser sempre analisado no final de cada iteração, calculando-se uma aproximação à solução com a precisão pretendida. 5 6 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR 1. A recolha de energia solar através da focagem de um campo plano de espelhos numa central de recolha foi estudada por Vant-Hull (1976). A equação para a concentração geométrica do factor C é dada por: C = π(h/cos(A))2F 0.5πD2(1 + sen(A) − 0.5cos(A)) em que A é o ângulo do campo, F é a cobertura da fracção do campo com espelhos, D é o diâmetro do colector e h é o comprimento do colector. Considerando h = 300, F = 0.8 e D = 14, calcule o ângulo positivo A inferior a π 25 para o qual a concentração do factor C é 1200. Utilize o método iterativo mais adequado e considere no critério de paragem ε1 = ε2 = 10−3 ou no máximo 3 iterações. Resolução: Mudança de variável: A → x e C → f e colocação da calculadora em modo radianos. Uma vez que se trata de uma função complicada, optou-se pelo Método da Secante, que não utiliza derivadas. São necessários dois pontos iniciais, x1 = 0 e x2 = π 25 ≈ 0.125664, uma vez que é dito que o ângulo A é positivo e inferior a π 25. Coloca-se a função na forma f(x) = 0, i.e.: f(x) = π(300/cos(x))20.8 0.5π142(1 + sen(x) − 0.5cos(x)) − 1200 = 0. Realiza-se a primeira iteração do método, através de (A.3), com k = 2: x3 = x2 − (x2 − x1)f(x2) f(x2) − f(x1) = 0.125664 − (0.125664 − 0) × (−13.845171) −13.845171 − 269.387755 = 0.119521 Em seguida, utiliza-se (A.4) para verificar o critério de paragem. Para tal calcula-se a segunda condição em primeiro lugar, uma vez que, se esta for falsa não há necessidade de calcular a primeira condição (estimativa do erro relativo ≤ ε1) e fica efectuado um cálculo necessário para a eventual iteração seguinte. 9 Assim, obtiveram-se os valores de x4 = 0.054701, f(x4) = −0.210997, x5 = 0.054718 e f(x5) = 0.000042. A verificação do critério de paragem é mostrada de seguida: |x5 − x4| |x5| = |0.054718 − 0.054701| |0.054718| = 0.000304 < 0.05 (Verdadeiro!) |f(x5)| = 0.000042 < 0.05 (Verdadeiro!) A solução óptima é encontrada ao fim de 3 iterações. O valor da taxa de juro corresponde a i∗ ≈ 0.054718. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = 20000*(x1*pot(1+x1,6))/(pot(1+x1,6)-1)-4000 Valores Iniciais x1 = 0.050000 x2 = 0.150000 Tolerâncias do Critério de Paragem e1 = 5.000000e-02 e2 = 5.000000e-02 Número Máximo de Iterações NMAX = 4 Iteração 1 xn = 0.054437 fxn = -3.562855 Iteração 2 xn = 0.054701 fxn = -0.210997 Iteração 3 xn = 0.054718 fxn = 0.000042 Número de Iterações Realizadas = 3 Solução xn = 0.054718 fxn = 0.000042 10 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR 3. A figura representa um pêndulo suspenso num tecto dum sala. O pêndulo baloiça-se de acordo com a seguinte expressão d = 80+90cos(π/3 t), t ≥ 0 em que d (cm) representa a distância até à parede de referência e depende do número de segundos t desde que o pêndulo foi posto em movimento. Calcule o instante de tempo t para o qual o pêndulo toca na parede da sala. Utilize o método de Newton, use para aproximação inicial t_1 = 4 e considere ε_1 = ε_2 = 10^{-3} ou no máximo 4 iterações. Apresente uma estimativa do erro relativo. Resolução: Mudança de variável: t → x e d → f e colocação da calculadora em modo radianos. Quando o pêndulo toca na parede, d = 0. Coloca-se a função na forma f(x) = 0, i.e., f(x) = 80 + 90cos(π/3 x) = 0 A derivada da função é f'(x) = −30πsen(π/3 x). Realiza-se a primeira iteração, aplicando a expressão (A.2), com k = 1 e x_1 = 4: x_2 = x_1 − f(x_1)/f'(x_1) = 4 − 35.0000/81.6210 = 3.5712 Verifica-se o critério de paragem, através de (A.4): |f(x_2)| = 5.6264 ≤ 0.001 (Falso!). Como a segunda condição não é verificada, não é necessário calcular a primeira (|x_2 − x_1| / |x_2| ≤ ε_1), prosseguindo-se com as iterações, que só finalizam quando as duas condições do critério de paragem forem verificadas em simultâneo. 11 Tem-se assim x3 = 3.46652, x4 = 3.4545 e x5 = 3.4544. Verificação do critério de paragem: |x5 − x4| |x5| = |3.4544 − 3.4545| |3.4544| = 0.000003 < 0.003 (Verdadeiro!) |f(x5)| = 0.0000 < 0.003 (Verdadeiro!) A solução é encontrada ao fim de 4 iterações. O instante de tempo é t∗ ≈ 3.4544 segundos. Estimativa do erro relativo: |x5 − x4| |x5| = |3.4544 − 3.4545| |3.4544| = 0.000003 . MATLAB: M-file: function [F,d] = f2_3(x) F = [80+90*cos(x*(pi/3))]; if nargout>1 d=-30*pi*sin(x*(pi/3)); end Comandos: >> x0=[4] >> options=optimset(’Jacobian’,’on’,’maxIter’,5,’TolX’,1.0e-3,’TolFun’,1.0e-3) >> [xsol,fsol,exitflag,output]=fsolve(’f2_3’,x0,options) Optimization terminated successfully: Relative function value changing by less than OPTIONS.TolFun xsol = 3.4544 fsol = 5.7489e-007 12 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR 4. O volume v de um líquido num tanque esférico de raio r está relacionado com a profundidade h do líquido da seguinte forma: v = πh2(3r − h) 3 . a) Calcule, utilizando um método que não recorre ao cálculo de derivadas, a pro- fundidade h, num tanque de raio r = 1 para um volume de 0.5. Utilize para aproximação inicial o intervalo [0.25, 0.5] e considere ε1 = ε2 = 10−2 ou no máximo 3 iterações. b) Repita os cálculos, nas mesmas condições da alínea anterior, mas utilizando para aproximação inicial o intervalo [2.5, 3]. Comente os resultados e analise a viabi- lidade da solução encontrada. Resolução: Mudança de variável: h → x e v → f. a) Utiliza-se o Método da Secante por este não usar derivadas, sendo necessários dois pontos iniciais, x1 = 0.25 e x2 = 0.5. Pretende-se resolver f(x) = 0.5. Coloca-se a expressão na forma f(x) = 0, i.e., πx2(3 × 1 − x) 3 − 0.5=0. Aplica-se (A.3), com k = 2: x3 = x2 − (x2 − x1)f(x2) f(x2) − f(x1) = 0.5 − (0.5 − 0.25) × 0.1545 0.1545 − (−0.3200) = 0.4186 . Emprega-se (A.4) para verificar o critério de paragem. |f(x3)| = | − 0.0263| ≤ 0.01 (Falso!) 13 Não se verifica a segunda condição de (A.4) e, como tal, não é necessário o cálculo da primeira. O processo iterativo continua até as duas condições serem verificadas simultaneamente. Obtêm-se os valores de x3 = 0.418601, x4 = 0.430450 e x5 = 0.431128. Verifica-se o critério de paragem: |x5 − x4| |x5| = |0.431128 − 0.430450| |0.431128| = 0.001573 < 0.01 (Verdadeiro!) |f(x5)| = 0.000015 < 0.01 (Verdadeiro!) A solução óptima é atingida ao fim de 3 iterações, com h∗ ≈ 0.431128. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = (pi/3)*pot(x1,2)*(3-x1)-0.5 Valores Iniciais x1 = 0.250000 x2 = 0.500000 Tolerâncias do Critério de Paragem e1 = 1.000000e-02 e2 = 1.000000e-02 Número Máximo de Iterações NMAX = 4 Iteração 1 xn = 0.418601 fxn = -0.026320 Iteração 2 xn = 0.430450 fxn = -0.001425 Iteração 3 xn = 0.431128 fxn = 0.000015 Número de Iterações Realizadas = 3 14 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR Solução xn = 0.431128 fxn = 0.000015 b) Aplica-se o mesmo método da alínea anterior, variando apenas os pontos iniciais, sendo x1 = 2.5 e x2 = 3. Obtêm-se os valores de x3 = 2.923606 e x4 = 2.944140. O critério de paragem é satisfeito no final de 2 iterações: |x4 − x3| |x4| = |2.944140 − 2.923606| |2.944140| = 0.006974 < 0.01 (Verdadeiro!) |f(x4)| = 0.007048 < 0.01 (Verdadeiro!) A profundidade do reservatório é de h∗ ≈ 2.944140. Este resultado, possível ma- tematicamente, não faz sentido no contexto do problema, i.e., não tem significado físico, uma vez que a profundidade do líquido é superior ao diâmetro do reservatório. A solução não é viável. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = (pi/3)*x1*x1*(3*1-x1)-0.5 Valores Iniciais x1 = 2.500000 x2 = 3.000000 Tolerâncias do Critério de Paragem e1 = 1.000000e-02 e2 = 1.000000e-02 Número Máximo de Iterações NMAX = 3 Iteração 1 xn = 2.923606 fxn = 0.183798 Iteração 2 xn = 2.944140 fxn = 0.007048 Número de Iterações Realizadas = 2 Solução xn = 2.944140 fxn = 0.007048 15 5. A figura representa um vulcão em erupção. A relação entre a distância y (milhas) percorrida pela lava e o tempo t (horas) é dada por: y = 7 (2 − 0.9t). Existe uma aldeia no sopé da montanha a uma distância de y = 10. O gabinete de protecção civil advertiu os moradores da aldeia de que a lava chegaria às suas casas em menos de 6 horas. Calcule utilizando um método iterativo que recorre ao cálculo de derivadas o instante de tempo em que a lava do vulcão atinge a aldeia. Considere ε1 = ε2 = 10−3 ou no máximo 3 iterações. Utilize nos cálculos 4 casas decimais. Nota: (ax)′ = ax ln(a) Resolução: Mudança de variável: t → x e y → f. Utiliza-se o Método de Newton, porque recorre ao cálculo de derivadas. A lava chega à aldeia quando y = 10. Coloca-se a expressão na forma f(x) = 0, i.e., 7(2 − 0.9x) − 10 = 0. A derivada da função é f ′(x) = −7 × 0.9x × ln(0.9). A primeira iteração realiza-se aplicando (A.2) com k = 1 e x1 = 6 (a informação de que a lava chegaria em menos de 6h é útil para a selecção da estimativa inicial). x2 = x1 − f(x1) f ′(x1) = 6 − 0.2799 0.3920 = 5.9278 Verifica-se o critério de paragem, com o auxílio da expressão (A.4). 16 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR |f(x2)| = 0.2515 ≤ 0.001 (Falso!). Uma vez que a segunda condição não é cumprida, não se realiza a primeira, e continua-se o método iterativo até ambas as condições serem verificadas em simultâ- neo. Obtêm-se os valores de x3 = 5.3114 e x4 = 5.3114. Verifica-se o critério de paragem: |x4 − x3| |x4| = |5.3114 − 5.3114| |5.3114| = 0 < 0.001 (Verdadeiro!) |f(x4)| = | − 1.4528 × 10−5| < 0.01 (Verdadeiro!) A solução é atingida ao fim de 3 iterações. O instante de tempo em que a lava do vulcão atinge a aldeia é t∗ ≈ 5.3114 horas. MATLAB: M-file: function [F,d] = f2_5(x) F = [7*(2-0.9^(x))-10]; if nargout>1 d=-7*0.9^(x)*log(0.9); end Comandos: >> x0=[6] >> options=optimset(’Jacobian’,’on’,’maxIter’,3,’TolX’,1.0e-3,’TolFun’,1.0e-3) >> [xsol, fsol,exitflag,output]=fsolve(’f2_5’,x0,options) Optimization terminated successfully: Relative function value changing by less than OPTIONS.TolFun xsol = 5.3114 fsol = -1.4528e-005 exitflag = 1 output = ... iterations: 3 ... 17 6. Em engenharia ambiental, a seguinte equação pode ser usada para calcular o nível de concentração de oxigénio c num rio, em função da distância x, medida a partir do local de descarga de poluentes: c(x) = 10 − 20(e−0.2x − e−0.75x). Calcule, usando um método que recorre ao cálculo de derivadas, a distância para a qual o nível de oxigénio desce para o valor 5. Utilize para aproximação inicial o valor x1 = 1.0 e considere ε1 = ε2 = 10−2 ou no máximo 3 iterações. Utilize nos cálculos 4 casas decimais. Resolução: Mudança de variável: c → f. Pretende-se resolver f(x) = 5. Coloca-se a expressão na forma f(x) = 0, i.e., 10 − 20(e−0.2x − e−0.75x) − 5 = 0. Utiliza-se o Método de Newton por este recorrer ao uso de derivadas, e aplica-se (A.2) com k = 1 e x1 = 1.0. A derivada é f ′(x) = 4e−0.2x − 15e−0.75x. x2 = x1 − f(x1) f ′(x1) = 1.0 − −1.9273 −3.8106 = 0.4942 . Usa-se (A.4) para verificar o critério de paragem. |f(x2)| = 0.6878 ≤ 0.01 (Falso!) Como a primeira condição do critério de paragem é falsa, não se calcula a segunda, e o processo iterativo continua até ambas as condições serem verdadeiras em simultâneo. 18 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR Obtêm-se os valores de x3 = 0.5964 e x4 = 0.6023. Verifica-se o critério de paragem: |x4 − x3| |x4| = |0.6023 − 0.5964| |0.6023| = 0.009832 < 0.01 (Verdadeiro!) |f(x4)| = |0.0001| < 0.01 (Verdadeiro!) A solução é atingida ao fim de 3 iterações. A distância para a qual o nível de oxigénio desce para o valor 5 é x∗ ≈ 0.6023. MATLAB: M-file: function [F,d] = f2_6(x) F = [10-20*(exp(-0.2*x)-exp(-0.75*x))-5]; if nargout>1 d=[4*exp(-0.2*x)-15*exp(-0.75*x)]; end Opções: >>x0=[1] >> options=optimset(’Jacobian’,’on’,’maxIter’,4,’TolX’,1.0e-2,’TolFun’,1.0e-2) Comando: >> [xsol, fsol,exitflag,output]=fsolve(’f2_6’,x0,options) Optimization terminated successfully: Relative function value changing by less than OPTIONS.TolFun xsol = 0.6023 fsol = 1.1350e-004 exitflag = 1 output = ... iterations: 4 ... 19 7. A velocidade v de um paraquedista é dada por v = g m c (1 − e −  c t m   ) em que g = 980. Para um paraquedista com m = 75000gr., calcule o coeficiente c (em gr./s), para o qual o paraquedista atinge a velocidade v = 3600cm/s ao fim de t = 6s Considere o intervalo inicial [10000, 15000]gr./s. Faça apenas duas iterações do método da secante. Apresente uma estimativa do erro relativo. Resolução: Mudança de variável: c → x e v → f. O problema é definido por f(x) = 3600. Coloca-se a expressão na forma f(x) = 0, i.e., f(x) = 980 × 75000 x (1 − e −   x × 6 75000   ) − 3600 = 0 Utiliza-se (A.3), com k = 2 e usando como pontos iniciais x1 = 10000 e x2 = 15000. x3 = x2 − (x2 − x1)f(x2) f(x2) − f(x1) = 15000 − (15000 − 10000) × (−175.8516) −175.8516 − 447.4321 = 13589.3132 De acordo com o enunciado, realiza-se a segunda iteração (k = 3). x4 = x3−(x3 − x2)f(x3) f(x3) − f(x2) = 13589.3132−(13589.3132 − 15000) × (−15.0185) −15.0185 − (−175.8516) = 13457.5837. 20 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR A solução obtida ao fim de 2 iterações do Método da Secante é c∗ ≈ 13457.5837gr./s. Estimativa do erro relativo: |x4 − x∗| |x∗| ≈ |x4 − x3| |x4| = |13457.5837 − 13589.3132| |13457.5837| = 0.0098 A estimativa do erro relativo é aproximadamente 0.0098 (0.98%). CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = (980*75000/x1)*(1-exp(-6*x1/75000))-3600 Valores Iniciais x1 = 10000.000000 x2 = 15000.000000 Tolerâncias do Critério de Paragem e1 = 1.000000e-08 e2 = 1.000000e-08 Número Máximo de Iterações NMAX = 2 Iteração 1 xn = 13589.313152 fxn = -15.018541 Iteração 2 xn = 13457.583681 fxn = 0.563691 Número máximo de iterações atingido. O processo ainda não convergiu. 21 8. A velocidade ascendente, v, de um foguetão pode ser calculada pela seguinte expressão: v = u ln(m_0 / (m_0 − q t)) − g t em que u é a velocidade relativa a que o combustível é expelido, m_0 é a massa inicial do foguetão no instante t = 0, q é a taxa de consumo de combustível e g é a aceleração da gravidade. Considerando u = 2200 m/s, g = 9.8m/s², m_0 = 1.6 × 10^5 Kg e q = 2680 Kg/s, calcule o tempo para o qual o foguetão atinge a velocidade v = 1000 m/s, sabendo que esse instante está entre 20 e 30 s. Utilize o método que achar mais adequado, com ε_1 = 10^{−2} e ε_2 = 10^{−1} ou no máximo 3 iterações. Resolução: Mudança de variável: t → x e v → f. O exercício é definido por f(x) = 1000. Coloca-se a expressão na forma f(x) = 0, i.e., 2200ln(160000 / (160000 − 2680x)) − 9.8x − 1000 = 0. Para evitar derivar, opta-se pelo Método da Secante. Aplica-se (A.3) com k = 2 e x_1 = 20 e x_2 = 30 (A informação de que esse instante está entre 20 e 30 permite estabelecer os valores iniciais). x_3 = x_2 − (x_2 − x_1)f(x_2) / (f(x_2) − f(x_1)) = 30 − (30 − 20) × 241.9514 / 241.9514 − (−298.4699) = 25.5229 Para verificar o critério de paragem, usa-se (A.4). |f(x_3)| = | − 23.0573| ≤ 0.1 (Falso!) 22 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR Não se verificando uma das condições, a outra não é calculada, e o processo iterativo continua até que as duas condições se verifiquem simultaneamente. Obtêm-se os valores de x3 = 25.5229, x4 = 25.9124 e x5 = 25.9426. Verifica-se o critério de paragem: |x5 − x4| |x5| = |25.9426 − 25.9124| |25.9426| = 0.0012 < 0.01 (Verdadeiro!) |f(x5)| = |0.0121| < 0.1 (Verdadeiro!) A solução é atingida ao fim de 3 iterações. O foguetão atinge a velocidade v = 1000m/s no instante t∗ ≈ 25.9426s. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = 2200*ln(1.6e+5/(1.6e+5-2680*x1))-9.8*x1-1000 Valores Iniciais x1 = 20.000000 x2 = 30.000000 Tolerâncias do Critério de Paragem e1 = 1.000000e-02 e2 = 1.000000e-01 Número Máximo de Iterações NMAX = 4 Iteração 1 xn = 25.522911 fxn = -23.057251 ... Número de Iterações Realizadas = 3 Solução xn = 25.942612 fxn = 0.012113 23 9. Considere a seguinte equação: C = M r [1−(1+r)−n] em que C é o capital emprestado, M é a mensalidade, r é a taxa de juro por cada período (expressa como uma fracção) e n é o número de anos. Uma pessoa pode pagar uma mensalidade de 1250 euros. Se pretender contrair um empréstimo de 10000 euros a 10 anos, qual é a taxa que poderá suportar? Use um método iterativo que não recorre à derivada, fazendo duas iterações e apre- sentando uma estimativa do erro relativo cometido. O valor da taxa deve estar entre 0.01 e 0.05. Resolução: Mudança de variável: r → x e C → f. A expressão no enunciado diz que 10000 = f(x). Coloca-se a expressão na forma f(x) = 0, i.e., 1250 x [1 − (1 + x)−10] − 10000 = 0. Usa-se o Método da Secante por não utilizar derivadas, aplicando (A.3) com k = 2 e x1 = 0.01 e x2 = 0.05 (estes valores foram retirados da informação de que o valor da taxa deve estar entre 0.01 e 0.05): x3 = x2 − (x2 − x1)f(x2) f(x2) − f(x1) = 0.05 − (0.05 − 0.01) × (−347.8313) (−347.8313) − 1839.1307 = 0.0436. De acordo com enunciado, realiza-se a segunda iteração com k = 3, f(x3) = −40.6643: x4 = x3 − (x3 − x2)f(x3) f(x3) − f(x2) = 0.0436 − (0.0436 − 0.05) × (−40.6643) (−40.6643) − (−347.8313) = 0.0428. Ao fim de duas iterações, o valor da taxa de juro é r∗ ≈ 0.0428. 24 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR Calcula-se uma estimativa do erro relativo, recorrendo a (A.4): |x4 − x3| |x4| = |0.0428 − 0.0436| |0.0428| = 0.0187. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = (1250/x1)*(1-pot((1+x1),-10))-10000 Valores Iniciais x1 = 0.010000 x2 = 0.050000 Tolerâncias do Critério de Paragem e1 = 1.000000e-08 e2 = 1.000000e-08 Número Máximo de Iterações NMAX = 2 Iteração 1 xn = 0.043638 fxn = -42.535793 Iteração 2 xn = 0.042752 fxn = 1.150584 Número máximo de iterações atingido. O processo ainda não convergiu. 25 10. Pela aplicação do Princípio de Arquimedes para determinação do calado de embarca- ções, pretende determinar-se a profundidade h correspondente ao equilíbrio tal que γsVs = γlVl(h) com γs = 918.35 kg/m3 (densidade do sólido), Vs = 1700m3 (volume do sólido), γl = 1.025kg/m3 (densidade do líquido) e Vl(h) volume do líquido deslocado, conforme a figura. Utilize o método de Newton para calcular o valor de h, supondo Vl(h) = h(h − 40)2. Utilize para aproximação inicial h(1) = 140 e ε1 = ε2 = 10−4, ou no máximo 3 iterações. Resolução: Mudança de variável: h → x e Vl → f. Na primeira expressão fazem-se as substituições das constantes e retira-se o valor de f(x): γsVs = γlVl(h) ⇔ Vl(h) = γsVs γl = 918.35 × 1700 1.025 = 1523117.073 Reescreve-se a equação na forma f(x) = 0, i.e., x(x − 40)2 − 1523117.073 = 0. Como se trata do Método de Newton, calcula-se a derivada: f ′(x) = 3x2−160x+1600. Aplica-se (A.2), com k = 1 e x1 = 140: x2 = x1 − f(x1) f ′(x1) = 140 − −123117.073 38000 = 143.2399. 26 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR Para verificar o critério de paragem, usa-se (A.4): |f(x2)| = |3602.0998| ≤ 0.0001 (Falso!) Como a segunda condição do critério de paragem não é cumprida, não se calcula a primeira, e o processo iterativo continua. Ao fim de 3 iterações a aproximação ainda não verifica o critério de paragem. No entanto, é atingido o número máximo de iterações estipulado e o processo iterativo pára, com x3 = 143.1504, f(x3) = 3.9016, x4 = 143.1503 e f(x4) = 3.2964e − 006. A profundidade é de h∗ ≈ 143.1503. MATLAB: M-File: function [F,d] = f2_10(x) F=[x*(x-40)^2-1523117.073]; if nargout>1 d=[3*x^2-160*x+1600]; end Opções: >> x0=[140] >> options=optimset(’Jacobian’,’on’,’maxIter’,3,’TolX’,1.0e-4,’TolFun’,1.0e-4) Comandos: >> [xsol, fsol,exitflag,output]=fsolve(’f2_10’,x0,options) Maximum number of iterations exceeded; increase options.MaxIter xsol = 143.1503 fsol = 0.5945 exitflag = 0 (indica que o processo iterativo ainda não convergiu) output = ... iterations: 3 ... 27 11. Uma das soluções para os resíduos de material nuclear é colocá-los em barris especiais que serão mais tarde depositados no fundo do oceano. Se os recipientes permanecerem intactos, a contaminação do ambiente circundante é mínima. Resolvendo as equações de movimento para os barris à medida que eles descem na água, chega-se à seguinte relação entre a velocidade de impacto, v, e a profundidade da água, D: D = 1 / k²g [W(W − B) ln(1 + kv / W − B) − Wkv], em que W é o peso dos barris, B é a sua flutuabilidade, g é a constante gravitacional e k é o coeficiente de atrito. A flutuabilidade dos barris pode ser determinada através do seu volume, sendo igual a 470. O coeficiente de atrito é determinado experimentalmente e é dado por k = 0.08. A constante gravitacional é g = 32 e o peso dos barris W = 527. a) Determine a velocidade de impacto v usando o método da secante, quando os barris são lançados numa zona cuja profundidade é D = −300. Utilize como aproximações iniciais v_1 = 40 e v_2 = 45, e no critério de paragem ε_1 = 0.05, ε_2 = 0.05 ou no máximo 2 iterações. b) Através de experiências, mostrou-se que os barris se danificam se a velocidade de impacto com o fundo do oceano for superior a 40. Na situação da alínea anterior, haverá risco de contaminação? Resolução: Mudança de variável: v → x e D → f. 28 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR a) O exercício é dado pela expressão −300 = f(x). Escreve-se a expressão na forma \(f(x) = 0\), i.e., \(\frac{1}{0.08^2 \times 32}\left[ 527(527 − 470) \ln\left(1 + \frac{0.08x}{527 − 470}\right) − 527 \times 0.08x\right] + 300 = 0\) Usando o Método da Secante, aplica-se (A.3), com k = 2, x1 = 40 e x2 = 45. x3 = x2 - \frac{(x2 - x1)f(x2)}{f(x2) - f(x1)} = 45 - \frac{(45 − 40)19.2249}{19.2249 − 77.1619} = 46.6591. Aplica-se (A.4) para verificação do critério de paragem. \(|f(x3)| = | −1.4162| \leq 0.05\) Uma vez que a segunda condição não é cumprida, não se calcula a primeira, fazendo-se nova iteração. Obtém-se o valor de x4 = 46.545284. Verifica-se o critério de paragem: \(\frac{|x4 − x3|}{|x4|} = \frac{|46.545284 − 46.659121|}{|46.545284|} = 0.002446 < 0.05\ (Verdadeiro!)\) \(|f(x4)| =|0.022420| < 0.05 (Verdadeiro!)\) A solução é atingida ao fim de 2 iterações. A velocidade de impacto é v* \approx 46.545284. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = (1/(pot(0.08,2)*32))*(527*(527-470)*ln(1+((0.08*x1)/(527-470)))-527*0.08*x1)+300 Valores Iniciais x1 = 40.000000 x2 = 45.000000 Tolerâncias do Critério de Paragem e1 = 5.000000e-02 e2 = 5.000000e-02 Número Máximo de Iterações 29 NMAX = 3 Iteração 1 xn = 46.659121 fxn = -1.416242 Iteração 2 xn = 46.545284 fxn = 0.022420 Número de Iterações Realizadas = 2 Solução xn = 46.545284 fxn = 0.022420 b) Nas condições da alínea anterior, há risco de contaminação porque a velocidade de impacto é v ≈ 46.545284, que é superior a 40, logo os barris são danificados. 30 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR 12. A pressão máxima, P, em Kg/mm2 que um cabo metálico suporta é dada por P(d) = 25d2 + ln(d) em que d é o diâmetro em mm. Determine o valor do diâmetro necessário para suportar uma pressão de 1.5 × 10−4 Kg/mm2. Sabendo que esse diâmetro pertence ao intervalo [0.2, 0.3], use o método iterativo que não recorre à informação da derivada para calcular uma aproximação com erro relativo inferior a 10−3. Use 6 casas decimais nos cálculos. Resolução: Mudança de variável: d → x e P → f. O problema é formulado como: 1.5 × 10−4 = 25x2 + ln(x). Coloca-se a expressão na forma f(x) = 0, i.e., 25x2 + ln(x) − 1.5 × 10−4 = 0. Para não utilizar derivadas, recorre-se ao Método da Secante, aplicando (A.3), com k = 2, x1 = 0.2 e x2 = 0.3 (sabe-se que o diâmetro pertence ao intervalo [0.2,0.3]). x3 = x2 − (x2 − x1)f(x2) f(x2) − f(x1) = 0.3 − (0.3 − 0.2)1.045877 1.045877 − (−0.609588) = 0.236823 No enunciado, o critério de paragem apenas refere a análise do erro relativo:. |x3 − x2| |x3| = | − 0.266770| ≤ 0.001 (Falso!) . 31 O processo iterativo continua. Obtêm-se os valores de x4 = 0.239064 e x5 = 0.239212. Verifica-se o critério de paragem: |x5 − x4| |x5| = |0.239212 − 0.239064| |0.239212| = 0.000619 < 0.001 (Verdadeiro!) |f(x5)| = |0.000006| A solução é atingida ao fim de 3 iterações. O diâmetro necessário é d∗ ≈ 0.239212 mm. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = 25*pot(x1,2)+ln(x1)-0.00015 Valores Iniciais x1 = 0.200000 x2 = 0.300000 Tolerâncias do Critério de Paragem e1 = 1.000000e-03 e2 = 1.000000e-03 Número Máximo de Iterações NMAX = 4 Iteração 1 xn = 0.236823 fxn = -0.038468 Iteração 2 xn = 0.239064 fxn = -0.002384 Iteração 3 xn = 0.239212 fxn = 0.000006 Número de Iterações Realizadas = 3 Solução xn = 0.239212 fxn = 0.000006 32 CAPÍTULO 2. SOLUÇÃO DE UMA EQUAÇÃO NÃO LINEAR 13. A função a(x) = 2.02x5 − 1.28x4 + 3.06x3 − 2.92x2 − 5.66x + 6.08 é utilizada num estudo do comportamento mecânico dos materiais, representando a(x) o comprimento da fissura e x(> 0) uma fracção do número de ciclos de propagação. Pretende-se saber para que valores de x a velocidade de propagação é nula. Utilize um método que não recorre ao cálculo de derivadas, usando como critério de paragem ε1 = ε2 = 10−2 ou no máximo três iterações. Resolução: Neste problema, a(x) representa o comprimento de uma fissura, mas pretende-se obter o ponto para o qual a velocidade de propagação da mesma é nula, i.e., a′(x) = 0. Para tal, é necessário encontrar o zero da derivada da função fornecida no enunciado, que corresponde a: a′(x) = 10.10x4 − 5.12x3 + 9.18x2 − 5.84x − 5.66. Mudança de variável: a′(x) → f(x). Coloca-se a função na forma f(x) = 0, i.e., 10.10x4−5.12x3+9.18x2−5.84x−5.66 = 0. Utiliza-se o Método da Secante, por não ser necessário o cálculo de derivadas, aplicando- se (A.3), com k = 2. Para tal, são necessários dois pontos iniciais, correspondentes ao intervalo onde se espera que a solução se encontre. Pela observação da figura, verifica-se que existem dois zeros para a velocidade, um em [−0.6, −0.4] e outro em 33 [0.8, 1]. Como um número de ciclos negativo não faria sentido, utiliza-se apenas o zero do intervalo positivo, com x1 = 0.8 e x2 = 1. x3 = x2 − (x2 − x1)f(x2) f(x2) − f(x1) = 1 − (1 − 0.8)2.66 2.66 − (−2.94128) = 0.905022. Usa-se (A.4) para verificação do critério de paragem. |f(x3)| = | − 0.445866| ≤ 0.01 (Falso!) O critério de paragem não foi cumprido, e o método iterativo prossegue até que ambas as condições se verifiquem em simultâneo. Obtêm-se os valores de x4 = 0.918657 e x5 = 0.920524. Verifica-se o critério de paragem: |x5 − x4| |x5| = |0.920524 − 0.918657| |0.920524| =< 0.01 (Verdadeiro!) |f(x5)| = |0.001327| < 0.01 (Verdadeiro!) A solução é obtida ao fim de 3 iterações. O ponto para o qual a velocidade de propagação da fissura é nula é x∗ ≈ 0.920524. CONUM: SECANT Solução de uma Equação Não Linear - Método da Secante Função f(x1) = 10.1*pot(x1,4)-5.12*pot(x1,3)+9.18*pot(x1,2)-5.84*x1-5.66 Valores Iniciais x1 = 0.800000 x2 = 1.000000 Tolerâncias do Critério de Paragem e1 = 1.000000e-02 e2 = 1.000000e-02 Número Máximo de Iterações NMAX = 4 Iteração 1 xn = 0.905022 fxn = -0.445874 ... Número de Iterações Realizadas = 3 Solução xn = 0.920524 fxn = 0.001327 Capítulo 3 Sistemas de equações lineares O objectivo deste capítulo é resolver exercícios que envolvem sistemas de n equações lineares do tipo:                a11x1 + a12x2 + . . . + a1nxn = b1 a21x1 + a22x2 + . . . + a2nxn = b2 ... ... ... ... = ... an1x1 + an2x2 . . . + annxn = bn Em termos matriciais fica Ax = b com A =   a11 a12 . . . a1n a21 a22 . . . a2n ... ... ... ... an1 an2 . . . ann   , x =   x1 x2 ... xn   e b =   b1 b2 ... bn   em que A ∈ ℜn×n é a matriz dos coeficientes, x ∈ ℜn é a solução do sistema e b ∈ ℜn é o termo independente. Os métodos para a sua resolução costumam classificar-se em duas categorias: os méto- dos directos e os métodos iterativos. Vão utilizar-se a Eliminação de Gauss com Pivotagem Parcial (EGPP) como método directo, e o método de Gauss-Seidel (A.5) como método ite- rativo com (A.7) no critério de paragem. Para este último, também se resolvem exercícios analisando as respectivas condições suficientes de convergência. 35 36 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES 1. Num determinado circuito eléctrico, as correntes i1, i2 e i3 passam através das impe- dâncias Z1, Z2 e Z3 e são dadas por:          i1 + i2 + i3 = 0 Z1i1 − Z2i2 = e1 − e2 Z2i1 − Z3i3 = e2 − e3 Se Z1 = 10, Z2 = 8, Z3 = 3, e1 − e2 = 65 e e2 − e3 = 120 : a) Calcule os valores das correntes i1, i2 e i3 por um método directo e estável. b) Calcule o determinante da matriz. c) Calcule a matriz inversa. Resolução: a) Mudança de variável: i → x. Substituindo as constantes, obtém-se:          x1 + x2 + x3 = 0 10x1 − 8x2 = 65 8x1 − 3x3 = 120 A =      1 1 1 10 −8 0 8 0 −3      b =      0 65 120      Procede-se à troca de linhas (−→ 1,2) porque o elemento de maior módulo da primeira coluna deve colocar-se na primeira linha, na primeira etapa.      1 1 1 | 0 10 −8 0 | 65 8 0 −3 | 120      −→ 1,2      10 −8 0 | 65 1 1 1 | 0 8 0 −3 | 120      37 1a etapa: Elemento pivot1 (a11): 10 (elemento de maior módulo da primeira coluna) Cálculo dos multiplicadores: m21 = − a21 pivot1 = − 1 10 = −0.1; m31 = − a31 pivot1 = − 8 10 = −0.8 O multiplicador m21 vai multiplicar a linha pivot (linha 1) e adicionar à linha 2. O multiplicador m31 vai multiplicar a linha pivot (linha 1) e adicionar à linha 3. Exemplo: −8 × (−0.1) + 1 = 1.8 A matriz ampliada obtida no final da 1a etapa é:      10 −8 0 | 65 0 1.8 1 | −6.5 0 6.4 −3 | 68      2a etapa: Trocam-se novamente as linhas (−→ 2,3), de modo a que o elemento de maior módulo da segunda coluna (da segunda linha para baixo) fique na posição a22.      10 −8 0 | 65 0 1.8 1 | −6.5 0 6.4 −3 | 68      −→ 2,3      10 −8 0 | 65 0 6.4 −3 | 68 0 1.8 1 | −6.5      Elemento pivot2 (a22): 6.4 (elemento de maior módulo da segunda coluna, a partir da segunda linha). Cálculo do multiplicador: m32 = − a32 pivot2 = −1.8 6.4 = −0.281250 O multiplicador m32 vai multiplicar a linha pivot (linha 2) e adicionar à linha 3.      10 −8 0 | 65 0 6.4 −3 | 68 0 0 1.84375 | −25.625      38 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES Assim, obteve-se o seguinte sistema, agora triangular, que se resolve por substituição inversa, ou seja, do fim para o início. Por exemplo, x3 = \frac{-25.625}{1.84375} = -13.898305, a seguir x2 = 4.110169 e por fim x1 = 9.788136. \left\{ \begin{array}{l} 10x_1 - 8x_2 = 65 \\ 6.4x_2 - 3x_3 = 68 \\ 1.84375x_3 = -25.625 \end{array} \right. Logo, os valores das diferentes correntes correspondem a i1 = 9.788136, i2 = 4.110169 e i3 = -13.898305. MATLAB: Comandos: >> A=[1 1 1; 10 -8 0;8 0 -3] A = 1 1 1 10 -8 0 8 0 -3 >> b=[0;65;120] b = 0 65 120 >> A\b ans = 9.7881 4.1102 -13.8983 b) Cálculo do determinante da matriz: \text{det}(A) = \text{det}(U) \times (-1)^t = \prod_{i=1,...,n}(u_{ii}) \times (-1)^t \ (t \text{ é o número de trocas de linhas}). A matriz U é a matriz triangular superior obtida no processo de eliminação de Gauss: \begin{pmatrix} 10 & -8 & 0 \\ 0 & 6.4 & -3 \\ 0 & 0 & 1.84375 \end{pmatrix} 39 Assim, det(A) = u11 × u22 × u33 × (−1)2 = 10 × 6.4 × 1.84375 × (−1)2 = 118. MATLAB: Comandos: >> A=[1 1 1; 10 -8 0;8 0 -3] >> det(A) ans = 118 c) Cálculo da matriz inversa (A|I) =      1 1 1 | 1 0 0 10 −8 0 | 0 1 0 8 0 −3 | 0 0 1      Aplica-se EGPP ao conjunto, efectuando as mesmas operações necessárias para o cál- culo da matriz triangular superior U (alínea a)) a partir da matriz (A|I). O resultado é o seguinte.      10 −8 0 | 0 1 0 0 6.4 −3 | 0 −0.8 1 0 0 1.84375 | 1 0.125 −0.28125      Para calcular a primeira coluna de A−1, resolve-se o sistema que tem como termo independente a 1a coluna da matriz da direita:      10 −8 0 | 0 0 6.4 −3 | 0 0 0 1.84375 | 1      De onde se obtém x1, x2 e x3 por substituição inversa, ficando calculada a primeira coluna de A−1. 1.84375x3 = 1 ⇔ x3 = 0.542373 6.4x2 − 3x3 = 0 ⇔ x2 = 0.254237 10x1 − 8x2 = 0 ⇔ x1 = 0.203390 40 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES Do mesmo modo, obtém-se para a segunda coluna:      10 −8 0 | 1 0 6.4 −3 | −0.8 0 0 1.84375 | 0.125      ⇔          x1 = 0.025424 x2 = −0.093220 x3 = 0.067797 E para a terceira coluna,      10 −8 0 | 0 0 6.4 −3 | 1 0 0 1.84375 | −0.28125      ⇔          x1 = 0.067797 x2 = 0.084746 x3 = −0.152542 Reunindo as soluções dos três sistemas anteriores surge a matriz inversa, A−1: A−1 =      0.203390 0.025424 0.067797 0.254237 −0.093220 0.084746 0.542373 0.067797 −0.152542      MATLAB: Comandos: >> A=[1 1 1; 10 -8 0;8 0 -3] >> inv(A) ans = 0.2034 0.0254 0.0678 0.2542 -0.0932 0.0847 0.5424 0.0678 -0.1525 41 2. Uma fábrica de tintas pretende utilizar as sobras de tinta de 4 tipos diferentes de tonalidades de tinta verde para criar uma tonalidade de verde mais popular. Uma unidade de medida (u.m.) da nova tinta será composta por x1 u.m. de tinta tipo 1, x2 u.m. de tinta tipo 2, x3 u.m. de tinta tipo 3 e x4 u.m. de tinta tipo 4. Cada u.m. de tinta nova é composta por 4 pigmentos que estão relacionados pelo seguinte sistema de equações lineares:                80x1+ 30x3+ 10x4 = 40 80x2+ 10x3+ 10x4 = 27 16x1+ 20x2+ 60x3+ 72x4 = 31 4x1+ 8x4 = 2 Os coeficientes da matriz representam a percentagem de pigmento em cada uma das 4 diferentes tonalidades de tinta verde, por exemplo, a tinta com a nova tonalidade deverá conter 31% de pigmento 3, sabendo que a tinta tipo 1 contem 16%, a tinta tipo 2 20%, a tinta tipo 3 60% e a tinta tipo 4 contem 72% do mesmo pigmento. a) Analisando apenas as condições suficientes de convergência, verifique se o método de Gauss-Seidel converge, quando aplicado a este sistema. b) Resolva o sistema de equações usando o método iterativo de Gauss-Seidel, utili- zando para aproximação inicial o ponto (0.5, 0.2, 0.2, 0)T e utilizando para cri- tério de paragem ε = 0.25 ou nmax = 2. Resolução: a) A =         80 0 30 10 0 80 10 10 16 20 60 72 4 0 0 8         b =         40 27 31 2         42 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES Análise das Condições Suficientes de Convergência: i) A matriz A é estrita e diagonalmente dominante? Não, porque |60| < |16| + |20| + |72| = 108. Nada se pode concluir quanto à convergência. ii) A matriz A é simétrica e definida positiva? A^T = \begin{pmatrix} 80 & 0 & 16 & 4 \\ 0 & 80 & 20 & 0 \\ 30 & 10 & 60 & 0 \\ 10 & 10 & 72 & 8 \end{pmatrix} \neq A ∴ A matriz A não é simétrica, pelo que não é necessário verificar se é definida positiva. Nada se pode concluir. iii) ||C_{GS}||_{1,\infty} < 1? (Cálculo da matriz de iteração C_{GS} (A.5)) D = \begin{pmatrix} 80 & 0 & 0 & 0 \\ 0 & 80 & 0 & 0 \\ 0 & 0 & 60 & 0 \\ 0 & 0 & 0 & 8 \end{pmatrix} L = \begin{pmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ -16 & -20 & 0 & 0 \\ -4 & 0 & 0 & 0 \end{pmatrix} U = \begin{pmatrix} 0 & 0 & -30 & -10 \\ 0 & 0 & -10 & -10 \\ 0 & 0 & 0 & -72 \\ 0 & 0 & 0 & 0 \end{pmatrix} D - L = \begin{pmatrix} 80 & 0 & 0 & 0 \\ 0 & 80 & 0 & 0 \\ 16 & 20 & 60 & 0 \\ 4 & 0 & 0 & 8 \end{pmatrix} (D - L)^{-1} = \begin{pmatrix} 0.0125 & 0 & 0 & 0 \\ 0 & 0.0125 & 0 & 0 \\ -0.0033 & -0.0042 & 0.0167 & 0 \\ -0.0063 & 0 & 0 & 0.1250 \end{pmatrix} C_{GS} = (D - L)^{-1}U = \begin{pmatrix} 0 & 0 & -0.3750 & -0.1250 \\ 0 & 0 & -0.1250 & -0.1250 \\ 0 & 0 & 0.1417 & -1.1250 \\ 0 & 0 & 0.1875 & 0.0625 \end{pmatrix} Cálculo da norma 1 da matriz C_{GS}: ||C_{GS}||_{1} = \max{\left(\sum_{i=1}^{n} |c_{i1}|, \sum_{i=1}^{n} |c_{i2}|, ..., \sum_{i=1}^{n} |c_{in}|\right)} = \max(0 + 0 + 0 + 0; 0 + 0 + 0 + 0; | - 0.3750| +| - 0.1250| + 0.1417 + 0.1875;| - 0.1250| +| - 0.1250| +| - 1.1250| + 0.0625) = \max(0; 0; 0.8292; 1.4375) = 1.4375. 43 ||C_{GS}||_{\infty} < 1 é falso, nada se pode concluir quanto à convergência. Cálculo da norma ∞ de C_{GS}: ||C_{GS}||_{\infty} = \max{\left(\sum_{j=1}^{n} |c_{1j}|, \sum_{j=1}^{n} |c_{2j}|, ..., \sum_{j=1}^{n} |c_{nj}|\right)} = \max(0+0+| -0.3750|+| -0.1250|; 0+ 0+| -0.1250| +| -0.1250|; 0 + 0 + 0.1417 + | - 1.1250|; 0 + 0 + 0.1875 + 0.0625) = \max(0.5; 0.25; 1.2667; 0.25) = 1.2667 ||C_{GS}||_{\infty} < 1 é falso, nada se pode concluir quanto à convergência. ||C_{GS}||_{1} = 1.4375 e ||C_{GS}||_{\infty} = 1.2667. Como nenhum destes valores é inferior a 1, nada se pode concluir. A análise das condições suficientes de convergência nada permite concluir acerca da convergência do Método de Gauss-Seidel. MATLAB: >> A=[80 0 30 10;0 80 10 10;16 20 60 72;4 0 0 8] >> D=diag(diag(A)) D = 80 0 0 0 0 80 0 0 0 0 60 0 0 0 0 8 >> L=tril(-A,-1) L = 0 0 0 0 0 0 0 0 -16 -20 0 0 -4 0 0 0 >> U=triu(-A,1) U = 0 0 -30 -10 0 0 -10 -10 0 0 0 -72 0 0 0 0 >> D-L 44 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES ans = 80 0 0 0 0 80 0 0 16 20 60 0 4 0 0 8 >> inv(D-L)*U ans = 0 0 -0.3750 -0.1250 0 0 -0.1250 -0.1250 0 0 0.1417 -1.1250 0 0 0.1875 0.0625 >> norm(inv(D-L)*U,1) ans = 1.4375 >> norm(inv(D-L)*U,inf) ans = 1.2667 b) x(1) =         0.5 0.2 0.2 0         Equação iterativa do Método de Gauss-Seidel (A.5): (D − L)xk+1 = Uxk + b 1a iteração (k = 1): (D − L)x(2) = Ux(1) + b         80 0 0 0 0 80 0 0 16 20 60 0 4 0 0 8         x(2) =         0 0 −30 −10 0 0 −10 −10 0 0 0 −72 0 0 0 0                 0.5 0.2 0.2 0         +         40 27 31 2         Nota: A matriz já é triangular (inferior), logo a resolução do sistema é por substituição 45 (80 & 0 & 0 & 0) \begin{pmatrix} x^{(2)} \end{pmatrix} = \begin{pmatrix} 34 \\ 25 \\ 31 \\ 2 \end{pmatrix} \Leftrightarrow \begin{pmatrix} x^{(2)} \end{pmatrix} = \begin{pmatrix} 0.425 \\ 0.3125 \\ 0.2992 \\ 0.0375 \end{pmatrix} x^{(2)} é obtido pela resolução do sistema por substituição directa. Por exemplo, a sua primeira componente é x_1 = \frac{34}{80} = 0.425. Verificação do critério de paragem (A.7): \begin{pmatrix} x^{(2)} \end{pmatrix} - \begin{pmatrix} x^{(1)} \end{pmatrix} = \begin{pmatrix} 0.425 \\ 0.3125 \\ 0.2992 \\ 0.0375 \end{pmatrix} - \begin{pmatrix} 0.5 \\ 0.2 \\ 0.2 \\ 0 \end{pmatrix} = \begin{pmatrix} -0.075 \\ 0.1125 \\ 0.0992 \\ 0.0375 \end{pmatrix} ||\begin{pmatrix} x^{(2)} \end{pmatrix} - \begin{pmatrix} x^{(1)} \end{pmatrix}|| = \sqrt{(-0.075)^2 + 0.1125^2 + 0.0992^2 + 0.0375^2 } = 0.1718 ||\begin{pmatrix} x^{(2)} \end{pmatrix}|| = \sqrt{0.425^2 + 0.3125^2 + 0.2992^2 + 0.0375^2} = 0.6076 = 0.2828 < 0.25 \text{ (Falso!) } O critério de paragem não é verificado, pelo que se continua para a segunda iteração. 2^a iteração (k=2): \begin{pmatrix} 80 & 0 & 0 & 0 \\ 0 & 80 & 0 & 0 \\ 16 & 20 & 60 & 0 \\ 4 & 0 & 0 & 8 \end{pmatrix} \begin{pmatrix} x^{(3)} \end{pmatrix} = \begin{pmatrix} 0 & 0 & -30 & -10 \\ 0 & 0 & -10 & -10 \\ 0 & 0 & 0 & -72 \\ 0 & 0 & 0 & 0 \end{pmatrix} \begin{pmatrix} 0.425 \\ 0.3125 \\ 0.2992 \\ 0.0375 \end{pmatrix} + \begin{pmatrix} 40 \\ 27 \\ 31 \\ 2 \end{pmatrix} \begin{pmatrix} 80 & 0 & 0 & 0 \\ 0 & 80 & 0 & 0 \\ 16 & 20 & 60 & 0 \\ 4 & 0 & 0 & 8 \end{pmatrix} \begin{pmatrix} x^{(3)} \end{pmatrix} = \begin{pmatrix} 30.649 \\ 23.633 \\ 28.3 \\ 2 \end{pmatrix} \equiv \begin{pmatrix} x^{(3)} \end{pmatrix} = \begin{pmatrix} 0.3831 \\ 0.2954 \\ 0.2710 \\ 0.0584 \end{pmatrix} Verificação do critério de paragem (A.7): 46 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES ||x(3) − x(2)|| ||x(3)|| = 0.0573 0.5576 = 0.1027 < 0.25 Verdadeiro! A estimativa do erro relativo é inferior a 0.25, e o processo iterativo termina. A aproximação à solução encontrada é: x∗ 1 ≈ 0.3831, x∗ 2 ≈ 0.2954, x∗ 3 ≈ 0.2710 e x∗ 4 ≈ 0.0584 u.m. de cada um dos tipos de tinta. CONUM: GAUSEI Sistemas de Equações Lineares - Equações de Gauss-Seidel Matriz A 80.000000 0.000000 30.000000 10.000000 0.000000 80.000000 10.000000 10.000000 16.000000 20.000000 60.000000 72.000000 4.000000 0.000000 0.000000 8.000000 Vector B 40.000000 27.000000 31.000000 2.000000 Vector Inicial 0.500000 0.200000 0.200000 0.000000 Tolerância do Critério de Paragem e = 2.500000e-01 Número Máximo de Iterações NMAX = 3 Matriz de Iteração do Método Gauss-Seidel 0.000000 0.000000 -0.375000 -0.125000 0.000000 0.000000 -0.125000 -0.125000 0.000000 0.000000 0.141667 -1.125000 0.000000 0.000000 0.187500 0.062500 O processo iterativo pode não convergir. Iteração 1 47 x[1] = 0.425000 x[2] = 0.312500 x[3] = 0.299167 x[4] = 0.037500 Iteração 2 x[1] = 0.383125 x[2] = 0.295417 x[3] = 0.271028 x[4] = 0.058437 Número de Iterações Realizadas = 2 Solução do Sistema 0.383125 0.295417 0.271028 0.058437 48 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES 3. Um engenheiro de Produção supervisiona a produção de quatro tipos de computa- dores. Existem quatro espécies de recursos necessários à produção: mão-de-obra, metais, plásticos e componentes electrónicos. As quantidades destes recursos, neces- sárias para produzir cada computador são: Mão de obra (h/comp.) Metais (Kg/comp.) Plásticos (Kg/comp.) Componentes (unid./comp.) 1 3 20 10 10 2 4 25 15 8 3 7 40 20 10 4 20 50 22 15 Considere um consumo diário de 504 h de mão-de-obra, 1970 Kg de metais, 970 Kg de plásticos e 601 componentes. a) Use um método directo e estável para calcular o número de computadores (número inteiro) de cada tipo produzidos por dia. b) Use o método iterativo de Gauss-Seidel, tomando como aproximação inicial x(1) = (9, 10, 12, 10). Apresente apenas os cálculos relativos às duas primeiras iterações, indicando uma estimativa do erro relativo. c) Comente os resultados obtidos, analisando as condições suficientes de convergên- cia. Resolução: a) Formulação do sistema de equações lineares:                3x1 + 4x2 + 7x3 + 20x4 = 504 20x1 + 25x2 + 40x3 + 50x4 = 1970 10x1 + 15x2 + 20x3 + 22x4 = 970 10x1 + 8x2 + 10x3 + 15x4 = 601 49 A =         3 4 7 20 20 25 40 50 10 15 20 22 10 8 10 15         b =         504 1970 970 601         Resolve-se o sistema linear por EGPP:         3 4 7 20 | 504 20 25 40 50 | 1970 10 15 20 22 | 970 10 8 10 15 | 601         −→ 1,2         20 25 40 50 | 1970 3 4 7 20 | 504 10 15 20 22 | 970 10 8 10 15 | 601         m21 = − 3 20 = −0.15, m31 = −10 20 = −0.5 e m41 = −10 20 = −0.5         20 25 40 50 | 1970 0 0.25 1 12.5 | 208.5 0 2.5 0 −3 | −15 0 −4.5 −10 −10 | −384         −→ 2,4         20 25 40 50 | 1970 0 -4.5 −10 −10 | −384 0 2.5 0 −3 | −15 0 0.25 1 12.5 | 208.5         m32 = − 2.5 −4.5 = 0.555556 e m42 = − 0.25 −4.5 = 0.055556         20 25 40 50 | 1970 0 −4.5 −10 −10 | −384 0 0 -5.555556 −8.555556 | −228.333333 0 0 0.444444 11.9444444 | 187.166667         m43 = − 0.444444 −5.555556 = 0.080000         20 25 40 50 | 1970 0 −4.5 −10 −10 | −384 0 0 −5.555556 −8.555556 | −228.333333 0 0 0 11.260000 | 169.900000         50 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES Resolve-se o sistema por substituição inversa, obtendo-se x∗ 1 = 10, x∗ 2 = 12, x∗ 3 = 18 e x∗ 4 = 15 computadores de cada tipo produzidos por dia. b) D =         3 0 0 0 0 25 0 0 0 0 20 0 0 0 0 15         L =         0 0 0 0 −20 0 0 0 −10 −15 0 0 −10 −8 −10 0         U =         0 −4 −7 −20 0 0 −40 −50 0 0 0 −22 0 0 0 0         D−L =         3 0 0 0 20 25 0 0 10 15 20 0 10 8 10 15         (D−L)−1 =         0.3333 0.0000 0.0000 −0.0000 −0.2667 0.0400 0 0.0000 0.0333 −0.0300 0.0500 −0.0000 −0.1022 −0.0013 −0.0333 0.0667         CGS = (D − L)−1U =         0 −1.3333 −2.3333 −6.6667 0 1.0667 0.2667 3.3333 0 −0.1333 0.9667 −0.2667 0 0.4089 0.7689 2.8444         Equação iterativa de Gauss-Seidel (A.5): (D − L)xk+1 = Uxk + b x(1) =         9 10 12 10         1a iteração(k = 1):         3 0 0 0 20 25 0 0 10 15 20 0 10 8 10 15         x(2) =         0 −4 −7 −20 0 0 −40 −50 0 0 0 −22 0 0 0 0                 9 10 12 10         +         504 1970 970 601         52 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES Nada se pode concluir quanto à convergência. ii) A é simétrica e definida positiva? AT =         3 20 10 10 4 25 15 8 7 40 20 10 20 50 22 15         ̸= A ∴ A não é simétrica. Nada se pode concluir quanto à convergência. iii) Análise de CGS CGS já foi calculado na alínea a) (CGS = (D − L)−1U). ||CGS||∞ < 1 ∨ ||CGS||1 < 1? ||CGS||∞ = 10.333323 < 1 (Falso!) ||CGS||1 = 13.111096 < 1 (Falso!) Nada se pode concluir quanto à convergência. A análise das condições suficientes de convergência não é conclusiva relativamente à convergência do Método de Gauss-Seidel para o sistema apresentado. 53 4. Considere o seguinte sistema de equações para determinar as concentrações c1, c2 e c3 (g/m3) numa série de 3 reactores como função da quantidade de massa à entrada de cada reactor (termo independente do sistema em g):          17c1 −2c2 −3c3 = 500 −5c1 +21c2 −2c3 = 200 −5c1 −5c2 +22c3 = 30 a) Analise as condições suficientes de convergência do método de Gauss-Seidel quando aplicado ao sistema. b) Aplique o método de Gauss-Seidel ao sistema, considerando como aproximação inicial o ponto (34, 19, 13) e ε1 = 0.0025 ou no máximo 2 iterações. Resolução: a) Análise das Condições Suficientes de Convergência: A =      17 −2 −3 −5 21 −2 −5 −5 22      b =      500 200 30      i) A matriz A é estrita e diagonalmente dominante? Sim, porque aii > Σn j=1,i̸=j|aij| i = 1, ..., n, ou seja |17| > | − 2| + | − 3|, |21| > | − 5| + | − 2| e |22| > | − 5| + | − 5|. Como uma das condições suficientes de convergência já é verificada, não é necessário verificar as outras, pelo que se conclui que o Método de Gauss-Seidel converge quando aplicado a este sistema. b) Equação iterativa do Método de Gauss-Seidel: (D − L)xk+1 = Uxk + b 54 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES D =      17 0 0 0 21 0 0 0 22      L =      0 0 0 5 0 0 5 5 0      U =      0 2 3 0 0 2 0 0 0      D − L =      17 0 0 −5 21 0 −5 −5 22      1a iteração (k = 1):      17 0 0 −5 21 0 −5 −5 22      x(2) =      0 2 3 0 0 2 0 0 0           34 19 13      +      500 200 30           17 0 0 −5 21 0 −5 −5 22      x(2) =      577 226 30      ⇔ x(2) =      33.941176 18.843137 13.360071      Verificação do Critério de Paragem (A.7): ||x(2) − x(1)|| ||x(2)|| = 0.397136 41.055557 = 0.009673 < 0.0025 (Falso!) O critério não é verificado e o processo iterativo prossegue. 2a iteração (k = 2):      17 0 0 −5 21 0 −5 −5 22      x(3) =      0 2 3 0 0 2 0 0 0           33.941176 18.843137 13.360071      +      500 200 30           17 0 0 −5 21 0 −5 −5 22      x(3) =      577.766487 226.720142 30      ⇔ x(3) =      33.986264 18.888165 13.380552      55 Verificação do Critério de Paragem: ||x(3) − x(2)|| ||x(3)|| = 0.066932 41.120167 = 0.001628 < 0.0025 (Verdadeiro!) O critério de paragem é cumprido. A solução é encontrada ao fim de 2 iterações e corresponde a c∗ 1 ≈ 33.986264, c∗ 2 ≈ 18.888165 e c∗ 3 ≈ 13.380552 g/m3. CONUM: GAUSEI Sistemas de Equações Lineares - Equações de Gauss-Seidel Matriz A 17.000000 -2.000000 -3.000000 -5.000000 21.000000 -2.000000 -5.000000 -5.000000 22.000000 Vector B Vector Inicial 500.000000 34.000000 200.000000 19.000000 30.000000 13.000000 Tolerância do Critério de Paragem Número Máximo de Iterações e = 2.500000e-03 NMAX = 3 Matriz de Iteração do Método Gauss-Seidel 0.000000 0.117647 0.176471 0.000000 0.028011 0.137255 0.000000 0.033104 0.071301 Iteração 1 x[1] = 33.941176 x[2] = 18.843137 x[3] = 13.360071 ... Número de Iterações Realizadas = 2 Solução do Sistema 33.986264 18.888165 13.380552 56 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES 5. Considere a figura representando um sistema de 4 molas ligadas em série sujeito a uma força F de 2000 Kg. Numa situação de equilíbrio, as equações força-balanço deduzidas definem inter- relações entre as molas:                k2(x2 − x1) = k1x1 k3(x3 − x2) = k2(x2−x1) k4(x4 − x3) = k3(x3−x2) F = k4(x4−x3) em que k1 = 150, k2 = 50, k3 = 75 e k4 = 225 são as constantes das molas (kg/s2). Analise as três condições suficientes de convergência do método de Gauss-Seidel e conclua sobre a convergência do método na resolução do sistema linear dado. Nota: a matriz de iteração CGS é CGS =   0.000 0.250 0.000 0.000 0.000 0.100 0.600 0.000 0.000 0.025 0.150 0.750 0.000 0.025 0.150 0.750   . Resolução: Substituição dos valores das constantes:                50(x2 − x1) = 150x1 75(x3 − x2) = 50(x2−x1) 225(x4 − x3) = 75(x3−x2) 2000 = 225(x4−x3) ⇔                −50x1 + 50x2 = 150x1 −75x2 + 75x3 = −50x1 + 50x2 −225x3 + 225x4 = −75x2 + 75x3 2000 = −225x3 + 225x4 57 Colocar em segundo membro os termos sem incógnitas:                −50x1 − 150x1 + 50x2 = 0 50x1 − 50x2 − 75x2 + 75x3 = 0 75x2 − 75x3 − 225x3 + 225x4 = 0 225x3 − 225x4 = 2000 ⇔                −200x1 + 50x2 = 0 50x1 − 125x2 + 75x3 = 0 75x2 − 300x3 + 225x4 = 0 225x3 − 225x4 = −2000 A =         −200 50 0 0 50 −125 75 0 0 75 −300 225 0 0 225 −225         b =         0 0 0 −2000         Análise das condições suficientes de convergência: i) A matriz A é estrita e diagonalmente dominante? Não, porque | − 125| > |50| + |75| ⇔ 125 > 125 (Falso!) ii) A matriz A é simétrica e definida positiva? AT =         −200 50 0 0 50 −125 75 0 0 75 −300 225 0 0 225 −225         = A ∴ A matriz A é simétrica. Verifica-se agora se a matriz A é definida positiva, calcu- lando os determinantes de todas as suas submatrizes. det(a11) = det(−200) > 0 (Falso!) Nada se pode concluir. iii) Como é fornecida a matriz CGS pode calcular-se imediatamente ||CGS||1 e ||CGS||∞ ||CGS||1,∞ < 1? ||CGS||1 = 1.5 e ||CGS||∞ = 0.925 Como uma das normas é menor que 1, conclui-se que o Método de Gauss-Seidel converge quando aplicado a este sistema. 58 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES 6. Uma equipa de três paraquedistas ligados por uma corda de peso desprezável é lan- çada em queda livre a uma velocidade v = 5 m/s conforme a figura. Considere os seguintes dados: Paraquedista Massa Coef. de resistência (i) (mi) (Kg) (ci) (Kg/s) 1 70 10 2 60 14 3 40 17 O sistema linear resultante permite calcular a tensão em cada secção da corda (R e T) e a aceleração da equipa (a).          m1g −c1v −T = m1a m2g −c2v +T −R = m2a m3g −c3v +R = m3a (considere g = 9.8 m/s2). O que poderia dizer acerca da convergência do método iterativo de Gauss-Seidel quando aplicado ao sistema? Justifique. Resolução: Mudança de variável: a → x1, R → x2 e T → x3 e substituição das constantes:          70 × 9.8 −10 × 5 −x3 = 70x1 60 × 9.8 −14 × 5 +x3 −x2 = 60x1 40 × 9.8 −17 × 5 +x2 = 40x1 59 Isolando os termos sem incógnitas no 2o membro:          70x1 +x3 = 636 60x1 +x2 −x3 = 518 40x1 −x2 = 307 A =      70 0 1 60 1 −1 40 −1 0      b =      636 518 307      Para determinar se o Método de Gauss-Seidel converge quando aplicado ao sistema, analisam-se as Condições Suficientes de Convergência: i) A matriz A é estrita e diagonalmente dominante? Não, porque |1| < |60| + | − 1|. ii) A matriz A é simétrica e definida positiva? AT =      70 60 40 0 1 −1 1 −1 0      ̸= A Logo, a matriz não é simétrica. iii) ||CGS||1,∞ < 1? CGS = (D − L)−1U D =      70 0 0 0 1 0 0 0 0      L =      0 0 0 −60 0 0 −40 1 0      U =      0 0 −1 0 0 1 0 0 0      D − L =      70 0 0 60 1 0 40 −1 0      60 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES D − L não tem inversa, o que é um dos pressupostos para a utilização do Método de Gauss-Seidel. Desta forma, conclui-se que este método não pode ser aplicado na resolução deste problema, na forma em que o sistema se encontra. No entanto, seria possível resolver o sistema através do método directo (EGPP). MATLAB: >> A=[70 0 1;60 1 -1;40 -1 0] >> D=diag(diag(A)) D = 70 0 0 0 1 0 0 0 0 >> L=tril(-A,-1) L = 0 0 0 -60 0 0 -40 1 0 >> U=triu(-A,1) U = 0 0 -1 0 0 1 0 0 0 >> D-L ans = 70 0 0 60 1 0 40 -1 0 >> inv(D-L) Warning: Matrix is singular to working precision. 61 7. Uma transportadora tem três tipos de camiões, Camião1, Camião2 e Camião3, que estão equipados para levar três tipos diferentes de máquinas de acordo com a seguinte tabela: máquina A máquina B máquina C Camião1 1 0 2 Camião2 1 1 1 Camião3 1 2 1 Por exemplo, o Camião1 pode levar uma máquina A, nenhuma máquina B e duas máquinas C. Supondo que cada camião vai com carga máxima, quantos camiões de cada tipo devemos enviar para transportar exactamente 12 máquinas A, 10 máquinas B e 16 máquinas C? a) Acha que consegue garantir a convergência do método iterativo de Gauss-Seidel na resolução deste problema? Justifique. b) Resolva o problema por um método directo e estável. Resolução: a) Camiãok → xk, com k = 1, 2, 3.          x1 + x2 + x3 = 12 x2 + 2x3 = 10 2x1 + x2 + x3 = 16 . A =      1 1 1 0 1 2 2 1 1      b =      12 10 16      Para verificar a convergência do método, analisam-se as Condições Suficientes de Convergência: 62 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES i) A é estrita e diagonalmente dominante? |1| > |1| + |1| (Falso!) A não é estrita e diagonalmente dominante. Nada se conclui quanto à convergência. ii) A é simétrica e positiva? AT =      1 0 2 1 1 1 1 2 1      ̸= A ∴ a matriz A não é simétrica, logo nada se conclui. iii) ||CGS||1,∞ < 1? CGS = (D − L)−1U D =      1 0 0 0 1 0 0 0 1      L =      0 0 0 0 0 0 −2 −1 0      U =      0 −1 −1 0 0 −2 0 0 0      D−L =      1 0 0 0 1 0 2 1 1      (D−L)−1 =      1 0 0 0 1 0 −2 −1 1      (D−L)−1U =      0 −1 −1 0 0 −2 0 2 4      ||CGS||1 = 7 e ||CGS||∞ = 6, logo nada se pode concluir quanto à convergência. Nenhuma das condições é satisfeita, pelo que nada se pode concluir sobre a conver- gência do Método de Gauss-Seidel quando aplicado à resolução deste sistema. b) Um método directo e estável é a Eliminação de Gauss com pivotagem parcial      1 1 1 | 12 0 1 2 | 10 2 1 1 | 16      −→ 1,3      2 1 1 | 16 0 1 2 | 10 1 1 1 | 12      63 m31 = −1 2 = −0.5      2 1 1 | 16 0 1 2 | 10 0 0.5 0.5 | 4      m32 = −0.5 1 = −0.5      2 1 1 | 16 0 1 2 | 10 0 0 −0.5 | −1      Resolvendo por substituição inversa obtém-se x1 = 4, x2 = 6 e x3 = 2. Devem utilizar-se 4 camiões do tipo 1, 6 camiões do tipo 2 e 2 camiões do tipo 3 para transportar a carga desejada. MATLAB: Comandos: >> A=[1 1 1;0 1 2;2 1 1] >> b=[12;10;16] >> A\b ans = 4 6 2 64 CAPÍTULO 3. SISTEMAS DE EQUAÇÕES LINEARES 8. Hoje é dia de S. Valentim. Dois rapazes pretendem comprar um ramo de flores, com rosas e tulipas, para oferecer às respectivas namoradas. Considere x1 o número de rosas e x2 o número de tulipas de cada ramo. O primeiro dos rapazes vai comprar o ramo à florista "Mil Pétalas" que cobra por cada rosa 2 euros e por cada tulipa 2 euros, gastando 10 euros. O segundo decide comprar o ramo na florista "Tudo em flor" mas esta ainda está inde- cisa relativamente ao preço a cobrar por cada rosa (consi- dere esse preço igual a "k") cobrando-lhe 3 euros por cada tulipa, ficando o ramo por 13 euros. a) Coloque o problema na forma de um sistema de equações lineares, em função de x1, x2 e k. b) Calcule justificando, o valor de k, por forma a garantir a convergência do mé- todo iterativo de Gauss-Seidel na resolução do sistema (use apenas as condições suficientes de convergência baseadas na matriz dos coeficientes, A). Resolução: a)    2x1 + 2x2 = 10 kx1 + 3x2 = 13 b) A =   2 2 k 3   b =   10 13   Analisam-se as condições suficientes de convergência: i) A matriz A é estrita e diagonalmente dominante? |2| > |2| (Falso!) 65 A matriz A não é estrita e diagonalmente dominante. Nada se conclui quanto à convergência. ii) A matriz A é simétrica e definida positiva? A matriz é simétrica se k = 2. Para ser definida positiva, é preciso verificar se todas as submatrizes de A têm determinante positivo: det(a11) = det(2) = 2 > 0 det   a11 a12 a21 a22   = det   2 2 2 3   = 2 > 0 Conclui-se que a matriz é definida positiva. Este último determinante deve ser cal- culado através de: det(A) = det(U) × (−1)t. Para k = 2, fica assegurada a convergência do método, pois verifica-se uma das condições suficientes de convergência. Capítulo 4 Sistemas de equações não lineares Neste capítulo vão ser resolvidos problemas envolvendo sistemas de n equações não lineares. Seja F : Ω −→ ℜn, em que Ω é um domínio em ℜn. O objectivo é determinar as soluções do sistema de equações F(x) = 0 ⇔          f1(x1, . . . , xn) = 0 . . . . . . fn(x1, . . . , xn) = 0 ou seja, o ponto x∗ = (x1, . . . , xn) que verifica simultaneamente todas as n equações. Neste tipo de sistemas é muito difícil demonstrar a existência e unicidade de zeros de F. Apenas em casos muito simples admite uma solução analítica pelo que, os métodos para a sua resolução são iterativos devido ao carácter não linear das suas equações. O método numérico utilizado é o de Newton (A.8) que envolve o cálculo da matriz do Jacobiano (A.9) e utiliza (A.10) como critério de paragem. 67 68 CAPÍTULO 4. SISTEMAS DE EQUAÇÕES NÃO LINEARES 1. Num colector solar, um balanço de energia na placa absorvente e na placa de vidro produz o seguinte sistema de equações não lineares nas temperaturas absolutas da placa absorvente (T1) e da placa de vidro (T2):    (T 4 1 + 0.06823T1) − (T 4 2 + 0.05848T2) = 0.01509 (T 4 1 + 0.05848T1) − (2T 4 2 + 0.11696T2) = 0 Usando como aproximação inicial à solução T (1) = (0.30, 0.30) faça uma iteração do método iterativo de Newton. Indique uma estimativa do erro relativo da aproximação. Resolução: Mudança de variável: T → x. O sistema é não linear nas variáveis. Colocam-se as funções na forma f(x) = 0:    f1(x1, x2) = 0 f2(x1, x2) = 0 ⇔    (x4 1 + 0.06823x1) − (x4 2 + 0.05848x2) − 0.01509 = 0 (x4 1 + 0.05848x1) − (2x4 2 + 0.11696x2) = 0 Constrói-se em seguida a matriz do Jacobiano.   ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2   =   4x3 1 + 0.06823 −4x3 2 − 0.05848 4x3 1 + 0.05848 −8x3 2 − 0.11696   Utiliza-se (A.8) com k = 1 e x(1) = (0.30, 0.30). f(0.30, 0.30) =   f1(0.30, 0.30) f2(0.30, 0.30)   =   −0.012165 −0.025644   J(0.30, 0.30) =   0.17623 −0.16648 0.16648 −0.33296   69 Resolve-se o seguinte sistema linear por EGPP, cujas incógnitas são ∆x1 e ∆x2:   0.17623 −0.16648 0.16648 −0.33296     ∆x1 ∆x2   (1) = −   −0.012165 −0.025644   ⇔ ⇔   0.17623 −0.16648 | 0.012165 0.16648 −0.33296 | 0.025644   ⇔   ∆x1 ∆x2   (1) =   −0.007065 −0.080551   A primeira iteração termina com o cálculo de x(2) = x(1) + ∆x(1).   x1 x2   (2) =   0.30 0.30   +   −0.007065 −0.080551   =   0.292935 0.219449   Para calcular uma estimativa do erro relativo, aplica-se (A.10): ||∆x(1)|| ||x(2)|| = 0.080860 0.366017 = 0.220919 MATLAB: M-file: function [F,d] = f4_1(x) F(1) = [(x(1)^4+0.06823*x(1))-(x(2)^4+0.05848*x(2))-0.01509]; F(2) = [(x(1)^4+0.05848*x(1))-(2*x(2)^4+0.11696*x(2))]; if nargout>1 d = [4*x(1)^3+0.06823 -4*x(2)^3-0.05848; 4*x(1)^3+0.05848 -8*x(2)^3-0.11696]; end Comandos: >> x=[0.30 0.30] >> options=optimset(’Jacobian’,’on’,’maxIter’,1) >> [xsol, fsol,exitflag,output]=fsolve(’f4_1’,x,options) Maximum number of iterations exceeded; increase options.MaxIter xsol =0.2929 0.2194 fsol = -0.0029 -0.0058 exitflag = 0 70 CAPÍTULO 4. SISTEMAS DE EQUAÇÕES NÃO LINEARES 2. Num jogo numa máquina de casino, as quantias em posse de dois jogadores são designadas por x1 e x2, e a quantia com que o casino fica é x3. Num determinado momento, essas quantias obedecem às equações fi = 0, i = 1, 2, 3, em que as fi são dadas por:          f1(x1, x2, x3) = x2 1 − x2 2 − 15 f2(x1, x2, x3) = 1 x2 1 + x2 2 − x3 f3(x1, x2, x3) = x1 + x2 + x3 − 500 Determine a quantidade com que cada um dos elementos fica nessa jogada. Considere a seguinte aproximação inicial (20, 20, 400) e apresente um resultado com erro relativo inferior a 0.01 ou no máximo 2 iterações. Resolução: Colocam-se as expressões na forma f(x) = 0, i.e.,          f1(x1, x2, x3) = 0 f2(x1, x2, x3) = 0 f3(x1, x2, x3) = 0 ⇔          x2 1 − x2 2 − 15 = 0 1 x2 1 + x2 2 − x3 = 0 x1 + x2 + x3 − 500 = 0 Constrói-se a matriz do Jacobiano,      ∂f1 ∂x1 ∂f1 ∂x2 ∂f1 ∂x3 ∂f2 ∂x1 ∂f2 ∂x2 ∂f2 ∂x3 ∂f3 ∂x1 ∂f3 ∂x2 ∂f3 ∂x3      =      2x1 −2x2 0 − 2 x3 1 2x2 −1 1 1 1      Utiliza-se (A.8) com k = 1 e x(1) = (20, 20, 400). f(20, 20, 400) =      f1(20, 20, 400) f2(20, 20, 400) f3(20, 20, 400)      =      −15 0.0025 −60      71 J(20, 20, 400) =      40 −40 0 −0.00025 40 −1 1 1 1      Resolve-se o sistema linear por EGPP:      40 −40 0 −0.00025 40 −1 1 1 1           ∆x1 ∆x2 ∆x3      (1) = −      −15 0.0025 −60      ⇔ ⇔      40 −40 0 | 15 −0.00025 40 −1 | −0.0025 1 1 1 | 60      ⇔      ∆x1 ∆x2 ∆x3      (1) =      1.794594 1.419594 56.785812      A iteração termina com o cálculo de x(2) = x(1) + ∆x(1).      x1 x2 x3      (2) =      20 20 400      +      1.794594 1.419594 56.785812      =      21.794594 21.419594 456.785812      Para estimar o erro relativo tem de se usar a primeira condição de (A.10). Calcula-se ||x(2) − x(1)||. ||x(2) − x(1)|| =      21.794594 21.419594 456.785812      −      20 20 400      =      1.794594 1.419594 56.785812      ||x(2) − x(1)|| ||x(2)|| = √ 1.7945942 + 1.4195942 + 56.7858122 √ 21.7945942 + 21.4195942 + 456.7858122 = 56.831895 457.806817 = 0.124139 Comparando a estimativa do erro relativo com 0.01, verifica-se que esta é superior e que a primeira condição do critério de paragem não é cumprida, pelo que se passa para a segunda iteração. 73 CONUM: NEWSIS Sistemas de Equações Não Lineares - Método de Newton Vector Funções pot(x1,2)-pot(x2,2)-15 1/(pot(x1,2))+pot(x2,2)-x3 x1+x2+x3-500 Matriz do Jacobiano 2*x1 -2*x2 0 -2/pot(x1,3) 2*x2 -1 1 1 1 Vector Inicial 20.000000 20.000000 400.000000 Tolerâncias do Critério de Paragem e1 = 1.000000e-02 e2 = 1.0 Número Máximo de Iterações NMAX = 2 ... Iteração 1 Iteração 2 x[1] = 21.794594 x[1] = 21.723359 x[2] = 21.419594 x[2] = 21.375248 x[3] = 456.785812 x[3] = 456.901392 fx[1] = 1.205321 fx[1] = 0.003108 fx[2] = 2.015301 fx[2] = 0.001967 fx[3] = 0.000000 fx[3] = 0.000000 Número de Iterações Realizadas = 2 Solução do Sistema Não Linear 21.723359 21.375248 456.901392 74 CAPÍTULO 4. SISTEMAS DE EQUAÇÕES NÃO LINEARES 3. A concentração de um poluente num lago depende do tempo t, e é dada por C(t) = 70eβt + 20eωt Efectuaram-se algumas medidas que foram registadas na seguinte tabela: t 1 2 C(t) 27.5702 17.6567 Utilize o método de Newton para determinar β e ω. Considere para aproximação inicial o ponto (β, ω)(1) = (−1.9, −0.15), efectue duas iterações e apresente uma estimativa do erro relativo. Resolução: Mudança de variável: β → x1 e ω → x2. Para construir as duas expressões que dão origem ao sistema, utiliza-se C(t) em cada um dos pontos, substituindo as respectivas constantes, obtendo-se:    C(1) = 27.5702 C(2) = 17.6567    70ex1 + 20ex2 − 27.5702 = 0 70e2x1 + 20e2x2 − 17.6567 = 0 ⇔    f1(x1, x2) = 0 f2(x1, x2) = 0 Constrói-se a matriz do Jacobiano.   ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2   =   70ex1 20ex2 140e2x1 40e2x2   Utiliza-se (A.8) com k = 1 e x(1) = (−1.9, −0.15). f(−1.9, −0.15) =   f1(−1.9, −0.15) f2(−1.9, −0.15)   =   0.113763 −1.274382   J(-1.9, -0.15) = \begin{pmatrix} 10.469803 & 17.214160 \\ 3.131908 & 29.632729 \end{pmatrix} Resolve-se o sistema linear por EGPP: \begin{pmatrix} 10.469803 & 17.214160 \\ 3.131908 & 29.632729 \end{pmatrix}\begin{pmatrix} \Delta x_1 \\ \Delta x_2 \end{pmatrix}^{(1)} = -\begin{pmatrix} 0.113763 \\ -1.274382 \end{pmatrix} \iff \begin{pmatrix} 10.469803 & 17.214160 & \ \vert -0.113763 \\ 3.131908 & 29.632729 & \vert & 1.274382 \end{pmatrix} \iff \begin{pmatrix} \Delta x_1 \\ \Delta x_2 \end{pmatrix}^{(1)} = \begin{pmatrix} -0.098732 \\ 0.053441 \end{pmatrix} Finalmente, calcula-se x^{(2)} = x^{(1)} + \Delta x^{(1)}. \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}^{(2)} = \begin{pmatrix} -1.9 \\ -0.15 \end{pmatrix} + \begin{pmatrix} -0.098732 \\ 0.053441 \end{pmatrix} = \begin{pmatrix} -1.998732 \\ -0.096559 \end{pmatrix} Estimativa do erro relativo (através da primeira condição de (A.10)): \frac{\|\Delta x^{(1)}\|}{\|x^{(2)}\|} = \frac{\sqrt{(-0.098732)^2 + 0.053441^2}}{\sqrt{(-1.998732)^2 + (-0.096559)^2}} = 0.056104 Conforme o enunciado, efectua-se a segunda iteração com x^{(2)} = (-1.998732, -0.096559): f(-1.998732, -0.096559) = \begin{pmatrix} f_1(-1.998732, -0.096559) \\ f_2(-1.998732, -0.096559) \end{pmatrix} = \begin{pmatrix} 0.074416 \\ 0.116342 \end{pmatrix} J(-1.998732, -0.096559) = \begin{pmatrix} 9.485490 & 18.159127 \\ 2.570700 & 32.975388 \end{pmatrix} Resolve-se o sistema linear por EGPP: \begin{pmatrix} 9.485490 & 18.159127 \\ 2.570700 & 32.975388 \end{pmatrix}\begin{pmatrix} \Delta x_1 \\ \Delta x_2 \end{pmatrix}^{(2)} = -\begin{pmatrix} 0.074416 \\ 0.116342 \end{pmatrix} \iff \begin{pmatrix} 9.485490 & 18.159127 & \vert -0.074416 \\ 2.570700 & 32.975388 & \vert -0.116342 \end{pmatrix} \iff \begin{pmatrix} \Delta x_1 \\ \Delta x_2 \end{pmatrix}^{(2)} = \begin{pmatrix} -0.001282 \\ -0.003428 \end{pmatrix} Finalmente, calcula-se x^{(3)} = x^{(2)} + \Delta x^{(2)}. \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}^{(3)} = \begin{pmatrix} -1.998732 \\ -0.096559 \end{pmatrix} + \begin{pmatrix} -0.001282 \\ -0.003428 \end{pmatrix} = \begin{pmatrix} -2.000014 \\ -0.099987 \end{pmatrix} Para a estimativa do erro relativo, usa-se a primeira condição de (A.10): \frac{\|\Delta x^{(2)}\|}{\|x^{(3)}\|} = \frac{\sqrt{(-0.001282)^2 + (-0.003428)^2}}{\sqrt{(-2.000014)^2 + (-0.099987)^2}} = 0.001828 Assim, obtêm-se os valores de \beta^* \approx -2.000014, \omega^* \approx -0.099987. CONUM: NEWSIS Sistemas de Equações Não Lineares - Método de Newton Vector Funções 70*exp(x1)+20*exp(x2)-27.5702 70*exp(2*x1)+20*exp(2*x2)-17.6567 Matriz do Jacobiano 70*exp(x1) 20*exp(x2) 140*exp(2*x1) 40*exp(2*x2) Vector Inicial -1.900000 -0.150000 Tolerâncias do Critério de Paragem e1 = 0.000000e+00 e2 = 0.000000e+00 Número Máximo de Iterações NMAX = 2 Iteração 1 Iteração 2 x[1] = -1.998732 x[1] = -2.000014 x[2] = -0.096559 x[2] = -0.099987 fx[1] = 0.074416 fx[1] = 0.000114 fx[2] = 0.116342 fx[2] = 0.000391 77 4. Para combater um vírus que infectou um grupo de indivíduos vai ser administrado um composto químico sintetizado com base em duas substâncias elementares x1 e x2. Sabe-se que se forem administrados α miligramas de composto a cada indivíduo, a concentração (mg/litro) de cada uma das substâncias elementares na circulação sanguínea é dada implicitamente (para α ∈ [0, 5]) pelo sistema de equações:    16x1 − cos(α(x2 − 2x1)) = 0 16x2 + 0.75 sen(α(−x2 − 3x1)) = 0 Para α = 1, determine x1 e x2 usando o método iterativo mais adequado. Use a seguinte aproximação inicial x(1) = (0.1, 0.01)T e termine o processo iterativo consi- derando ε2 = 0.05 (1 iteração). Resolução: Colocação da calculadora em modo Radianos. Substituindo no sistema α por 1 fica:    f1(x1, x2) = 0 f2(x1, x2) = 0 ⇔    16x1 − cos(x2 − 2x1) = 0 16x2 + 0.75 sen(−x2 − 3x1) = 0 Constrói-se a matriz do Jacobiano:   ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2   =   16 − 2sen(x2 − 2x1) sen(x2 − 2x1) −2.25cos(−x2 − 3x1) 16 − 0.75cos(−x2 − 3x1)   Utiliza-se (A.8), com k = 1 e (x1, x2)(1) = (0.1, 0.01). f(0.1, 0.01) =   f1(0.1, 0.01) f2(0.1, 0.01)   =   0.617996 −0.068794   J(0.1, 0.01) =   16.377718 −0.188859 −2.142751 15.285750   78 CAPÍTULO 4. SISTEMAS DE EQUAÇÕES NÃO LINEARES Resolve-se o sistema linear por EGPP:   16.377718 −0.188859 −2.142751 15.285750     ∆x1 ∆x2   (1) = −   0.617996 −0.068794   ⇔ ⇔   16.377718 −0.188859 | −0.617996 −2.142751 15.285750 | 0.068794   ⇔   ∆x1 ∆x2   (1) =   −0.037743 −0.000790   Calcula-se a nova aproximação, x(2) = x(1) + ∆x(1):   x1 x2   (2) =   0.062257 0.009210   Segue-se a verificação do critério de paragem, através de (A.10). f(x(2)) =   f1(0.062257, 0.009210) f2(0.062257, 0.009210)   =   0.002752 0.001313   ||f(x(2))|| = √ 0.0027522 + 0.0013132 = 0.003049 ≤ 0.05 (Verdadeiro!) O critério de paragem, conforme enunciado do exercício, envolve apenas a segunda condição de (A.10). O processo iterativo termina e a concentração das substâncias elementares é x∗ 1 ≈ 0.062087 e x∗ 2 ≈ 0.009099 (mg/litro). CONUM: NEWSIS Sistemas de Equações Não Lineares - Método de Newton Vector Funções 16*x1-cos(x2-2*x1) 16*x2+0.75*sen(-x2-3*x1) Matriz do Jacobiano 16-2*sen(x2-2*x1) sen(x2-2*x1) -2.25*cos(-x2-3*x1) 16-0.75*cos(-x2-3*x1) 79 Vector Inicial 0.100000 0.010000 Tolerâncias do Critério de Paragem e1 = 1.0 e2 = 5.000000e-02 Número Máximo de Iterações NMAX = 2 Iteração 1 x[1] = 0.062257 x[2] = 0.009210 fx[1] = 0.002751 fx[2] = 0.001309 Iteração 2 x[1] = 0.062087 x[2] = 0.009099 fx[1] = 0.000000 fx[2] = 0.000000 Número de Iterações Realizadas = 2 Solução do Sistema Não Linear 0.062087 0.009099 80 CAPÍTULO 4. SISTEMAS DE EQUAÇÕES NÃO LINEARES 5. Duas estações eléctricas vão fornecer energia a uma certa região da forma mais eco- nómica possível. O custo total de operação das duas estações é dado por: f(x1, x2) = 0.1+0.01x1x2+0.15x4 2+0.01x4 1−0.25(x1+x2−100) em que x1 é a energia fornecida pela primeira estação e x2 é a energia fornecida pela segunda estação. Determine os valores de x1 e x2 por forma a minimizar o custo total de operação das duas estações. Utilize como aproximação inicial o ponto (2.0, 0.5)T e ε1 = ε2 = 0.2 (uma iteração). Resolução: Neste exercício, pretende-se encontrar o mínimo do custo total de operação, que é obtido através dos zeros das primeiras derivadas da função f(x) em ordem a cada uma das variáveis. Realizando as derivadas parciais em ordem a x1 e x2, escreve-se o sistema na forma f(x) = 0, i.e.:    ∂f ∂x1 = 0 ∂f ∂x2 = 0 ⇔    f1(x1, x2) = 0 f2(x1, x2) = 0 ⇔    0.01x2 + 0.04x3 1 − 0.25 = 0 0.01x1 + 0.60x3 2 − 0.25 = 0 Constrói-se em seguida a matriz do Jacobiano do sistema anterior:   ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2   =   0.12x2 1 0.01 0.01 1.80x2 2   Utiliza-se (A.10) com k = 1 e x(1) = (2.0, 0.5). f(2.0, 0.5) =   f1(2.0, 0.5) f2(2.0, 0.5)   =   0.075 −0.155   J(2.0,0.5)= \begin{pmatrix} 0.48 & 0.01 \\ 0.01 & 0.45 \end{pmatrix} Resolve-se o sistema linear por EGPP: \begin{pmatrix} 0.48 & 0.01 \\ 0.01 & 0.45 \end{pmatrix}\begin{pmatrix} \Delta x_1 \\ \Delta x_2 \end{pmatrix}^{(1)} = -\begin{pmatrix} 0.075 \\ -0.155 \end{pmatrix} \begin{pmatrix} \Delta x_1 \\ \Delta x_2 \end{pmatrix}^{(1)} = \begin{pmatrix} -0.163502 \\ 0.348078 \end{pmatrix} A primeira iteração termina com o cálculo do novo valor de x, x^{(2)} = x^{(1)} + \Delta x^{(1)}. \begin{pmatrix} x_1 \\ x_2 \end{pmatrix}^{(2)} = \begin{pmatrix} 2.0 \\ 0.5 \end{pmatrix} + \begin{pmatrix} -0.163502 \\ 0.348078 \end{pmatrix} = \begin{pmatrix} 1.836498 \\ 0.848078 \end{pmatrix} Para verificar o critério de paragem, utiliza-se (A.10). f(x^{(2)}) = \begin{cases} f_1(1.836498, 0.848078) = 0.006241 \\ f_2(1.836498, 0.848078) = 0.134346 \end{cases} ||f(x^{(2)})|| = \sqrt{0.006241^2 + 0.134346^2} = 0.134491 \leq 0.2 \text{ (Verdadeiro!)} Verificando a outra condição: \frac{\|\Delta x^{(1)}\|}{\|x^{(2)}\|} = \frac{\sqrt{(-0.163502)^2 + 0.348078^2}}{\sqrt{1.836498^2 + 0.848078^2}} = 0.384566 = 0.190110 \leq 0.2 \text{ (Verdadeiro!)} O critério de paragem é verificado. O processo iterativo termina: x^*_1 \approx 1.836498 e x^*_2 \approx 0.848078. CONUM: NEWSIS Sistemas de Equações Não Lineares - Método de Newton Vector Funções 82 CAPÍTULO 4. SISTEMAS DE EQUAÇÕES NÃO LINEARES 0.01*x2+0.04*pot(x1,3)-0.25 0.01*x1+0.60*pot(x2,3)-0.25 Matriz do Jacobiano 0.12*pot(x1,2) 0.01 0.01 1.80*pot(x2,2) Vector Inicial 2.000000 0.500000 Tolerâncias do Critério de Paragem e1 = 2.000000e-01 e2 = 2.000000e-01 Número Máximo de Iterações NMAX = 1 Iteração 1 x[1] = 1.836498 x[2] = 0.848078 fx[1] = 0.006241 fx[2] = 0.134346 Número de Iterações Realizadas = 1 Solução do Sistema Não Linear 1.836498 0.848078 83 6. Considere a seguinte figura de uma viga em balanço: Um modelo de elementos finitos de uma viga em balanço sujeita a carga e momentos é obtido pela optimização de f(x, y) = 5x2 − 5xy + 2.5y2 − x − 1.5y, em que x e y são o deslocamento e o momento da extremidade, respectivamente. Calcule os valores de x e y que minimizam f(x, y), utilizando o método iterativo de Newton. Para aproximação inicial use (1, 1) e ε1 = ε2 = 10−6 ou no máximo duas iterações. Comente os resultados. Resolução: Mudança de variável: x → x1 e y → x2. Neste exercício, pretende-se encontrar o ponto que minimiza a função, que cor- responde aos zeros das primeiras derivadas em ordem às variáveis. Para tal, as funções igualadas a zero serão as derivadas parciais da função dada:    ∂f ∂x1 = f1(x1, x2) = 0 ∂f ∂x2 = f2(x1, x2) = 0 ⇔    10x1 − 5x2 − 1 = 0 −5x1 + 5x2 − 1.5 = 0 Constrói-se em seguida o Jacobiano:   ∂f1 ∂x1 ∂f1 ∂x2 ∂f2 ∂x1 ∂f2 ∂x2   =   10 −5 −5 5   Utiliza-se (A.8) com k = 1 e x(1) = (1, 1). f(1, 1) =   4 −1.5   84 CAPÍTULO 4. SISTEMAS DE EQUAÇÕES NÃO LINEARES J(1, 1) =   10 −5 −5 5   Resolve-se o sistema linear por EGPP:   10 −5 −5 5     ∆x1 ∆x2   (1) = −   4 −1.5   ⇔ ⇔   10 −5 | −4 −5 5 | 1.5   ⇔   ∆x1 ∆x2   (1) =   −0.5 −0.2   A primeira iteração termina com o cálculo do novo valor de x, x(2): x(2) = x(1) + ∆x(1) =   1 1   = +   −0.5 −0.2   =   0.5 0.8   O sistema que se resolveu é linear (reparar que a matriz do Jacobiano é constante) e portanto a solução exacta é obtida no final da primeira iteração. f(x(2)) =   0 0   e x∗ 1 = 0.5 , x∗ 2 = 0.8 . MATLAB: M-file: function [F,d] = f4_6(x) F(1) = [10*x(1)-5*x(2)-1]; F(2) = [-5*x(1)+5*x(2)-1.5]; if nargout>1 d = [10 -5;-5 5]; end 85 Comandos: >> x=[1 1] >> options=optimset(’Jacobian’,’on’,’TolX’, 1.0e-6,’TolFun’,1.0e-6,’maxIter’,2) >> [xsol, fsol,exitflag,output]=fsolve(’f4_6’,x,options) Optimization terminated successfully: First-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected xsol = 0.5000 0.8000 fsol = 1.0e-014 * (-0.1776 -0.1776) exitflag = 1 98 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA • (x_1 − x_0)M_0+2(x_2−x_0)M_1+(x_2−x_0)M_2 = 2(4− 1)M_1 + (4− 3)M_2 = \frac{6}{4−3}(2.25− 1.75)− \frac{6}{3−1}(1.75− 1) \Leftrightarrow 2M_0 • (x_2−x_1)M_1+2(x_3−x_1)M_2+(x_3−x_1)M_3 = 2(6−3)M_2+(6−4)M_3 = \frac{6}{6−4}(3.15−2.25)− \frac{6}{4−3}(2.25− 1.75) \Leftrightarrow 1M_1 Juntando as 4 equações anteriores na ordem dos nós x_0, x_1, x_2 e x_3, obtém-se o seguinte sistema linear nas incógnitas M_0, M_1, M_2 e M_3: \begin{pmatrix}4 & 2 & 0 & 0 | 0 & 0.75\\2 & 6 & 1 & 0 | 0 & 0.75\\0 & 1 & 6 & 2 | 0 &−0.3\\0 & 0 & 2 & 4|−1.8\end{pmatrix} Resolve-se o sistema linear por EGPP: M_0 = 0.160938, M_1 = 0.053125, M_2 = 0.109375 e M_3 =−0.504688 O ponto interpolador x = 1.5 localiza-se no 1° segmento [1,3]. Utiliza-se (A.14) com i = 1: s_3^1(x) = \frac{M_0}{6(x_1− x_0)}(x_1− x_0)^3 + \frac{M_1}{6(x_1− x_0)}(x− x_0)^3 + \left[\frac{f_0}{x_1− x_0}− \frac{M_0(x_1− x_0)}{6}\right](x− x_0) s_3^1(x) = 0.013412(3− x)^3 + 0.004427(x− 1)^3 + 0.446354(3− x) + 0.857292(x− 1) s_3^1(1.5) = 1.143996 \approx f(1.5) b) A aproximação à velocidade corresponde à derivada de s_3^1(x): s_3'^1(x) = −0.040236(3− x)^2 + 0.013281(x− 1)^2 s_3'^1(1.5) = 0.323727rad/s \approx Velocidade Capítulo 5 Interpolação Numérica A interpolação é utilizada principalmente em duas situações: quando a expressão de f(x) não é conhecida (apenas se conhecem os valores numéricos de f(x) para um conjunto discreto de pontos), ou quando se conhece a expressão da função f(x) mas ela é muito complicada, pelo que é útil aproximá-la por um modelo matemático mais simples do tipo polinomial. O problema geral da interpolação polinomial consiste em, dados n + 1 pontos distin- tos x0, x1, . . . , xn e respectivos valores y0, y1, . . . , yn em que yi = f(xi), determinar um polinómio pn(x) de grau n tal que pn(xi) = yi, i = 0, . . . , n (polinómio colocativo). Vão ser resolvidos exercícios aplicando interpolação simples utilizando o polinómio in- terpolador de Newton (A.12) também conhecido por diferenças divididas, estimando o respectivo erro (A.13). Quando se utilizam polinómios de grau muito elevado, estes assumem um compor- tamento muito oscilatório, não reflectindo o comportamento da função que pretendem estimar. A ideia subjacente à interpolação segmentada ou por splines (A.14, A.15) é a de utilizar vários polinómios de grau baixo em diferentes intervalos de pontos. São resolvidos exercícios utilizando splines cúbicas que, para além de possuírem derivadas contínuas até à segunda ordem, são fáceis de construir e apresentam um comportamento estável. Os exercícios contemplam a spline cúbica natural (A.16) e a spline cúbica completa (A.17). É também estimado o erro de truncatura (A.19). 87 88 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA 1. A tabela seguinte apresenta a população dos Estados Unidos da América (em mi- lhões) de 1940 a 1980. Ano 1940 1950 1960 1970 1980 População 132.165 151.326 179.323 203.302 226.542 a) Construa o polinómio interpolador de Newton de grau 4 para estimar a população no ano 1965. b) A população em 1930 foi 123.203 Qual a precisão do valor calculado na alínea a)? Resolução: Mudança de variável: ano → x e população → f(x). a) Se o polinómio é de grau 4 (n) são necessários 5 (n + 1) pontos. Usa-se (A.12). p4(x) = f0 + (x − x0)[x0, x1] + (x − x0)(x − x1)[x0, x1, x2] + (x − x0)(x − x1)(x − x2)[x0, x1, x2, x3] + (x − x0)(x − x1)(x − x2)(x − x3)[x0, x1, x2, x3, x4] Diferenças divididas: xi fi dd1 dd2 dd3 dd4 x0 = 1940 132.165 1.916100 x1 = 1950 151.326 0.044180 2.799700 −0.002142 x2 = 1960 179.323 −0.020090 0.000067 2.397900 0.000547 x3 = 1970 203.302 −0.003695 2.324000 x4 = 1980 226.542 Assim, p4(x) = 132.165+(x−1940)1.916100+(x−1940)(x−1950)0.044180+(x−1940)(x− 1950)(x − 1960) × (−0.002142) + (x − 1940)(x − 1950)(x − 1960)(x − 1970)0.000067 89 Para estimar a população no ano 1965, calcula-se p4(1965) ≈ f(1965): p4(1965) = 132.165 + (1965 − 1940)1.916100 + (1965 − 1940)(1965 − 1950)0.044180 + (1965 − 1940)(1965 − 1950)(1965 − 1960) × (−0.002142) + (1965 − 1940)(1965 − 1950)(1965 − 1960)(1965 − 1970)0.000067 = 191.987930 Em 1965, existiam aproximadamente 191.987930 milhões de habitantes nos Estados Unidos da América. b) Para calcular a precisão do valor calculado na alínea a), é necessário introduzir o valor de x = 1930 na tabela das diferenças divididas. O local de colocação do novo ponto é indiferente (a dd5 será sempre a mesma). Assim: xi fi dd1 dd2 dd3 dd4 dd5 1940 132.165 1.9161 1950 151.326 0.04418 2.7997 −0.002142 1960 179.323 −0.02009 0.000067 2.3979 0.000547 0.000002 1970 203.302 −0.003695 0.000044 2.324 −0.000338 1980 226.542 0.006431 2.06678 1930 123.203 Para calcular o erro utiliza-se a expressão (A.13), ou seja: |R4(x)| ≤ |(x − x0)(x − x1)(x − x2)(x − x3)(x − x4)||dd5| O erro de truncatura para o valor da alínea a), x = 1965, é de: |Rn(1965)| ≤ |(1965−1940)(1965−1950)(1965−1960)(1965−1970)(1965−1980)||0.000002| = 0.28125 90 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA CONUM: DIFDIV Interpolação Polinomial - Diferenças Divididas Pontos Valores da Função 1940.000000 132.165000 1950.000000 151.326000 1960.000000 179.323000 1970.000000 203.302000 1980.000000 226.542000 Ponto Interpolador = 1965.000000 Grau do Polinómio = 4 Interpolação Directa Tabela das Diferenças Diferenças de Ordem 1 1.916100 2.799700 2.397900 2.324000 Diferenças de Ordem 2 0.044180 -0.020090 -0.003695 Diferenças de Ordem 3 -0.002142 0.000547 Diferenças de Ordem 4 0.000067 Aproximação ao Valor da Função no Ponto Interpolador: 191.987930 91 2. A tabela seguinte apresenta a velocidade de queda de um paraquedista em função do tempo: tempo (s) 1 3 5 7 20 vel (cm/s) 800 2310 3090 3940 8000 a) Estime o valor da velocidade no instante de tempo t = 10s, utilizando um poli- nómio interpolador de grau 3. b) Calcule uma aproximação do erro cometido na alínea anterior. Resolução: Mudança de variável: tempo → x e vel → f(x). a) O polinómio interpolador de grau 3 necessita de 4 pontos e é dado por (A.12): p3(x) = f0 + (x − x0)[x0, x1] + (x − x0)(x − x1)[x0, x1, x2] + (x − x0)(x − x1)(x − x2)[x0, x1, x2, x3] Como seleccionar os 4 pontos a utilizar na construção do polinómio? O ponto inter- polador é x = 10 - são necessários os pontos imediatamente à esquerda e à direita, i.e., x = 7 e x = 20. Tendo já escolhido 2 pontos, falta ainda escolher os outros 2. Estes deverão estar o mais próximo possível do valor a interpolar (x = 3 e x = 5). xi fi dd1 dd2 dd3 x0 = 3 2310 390 x1 = 5 3090 8.75 425 −0.956637 x2 = 7 3940 −7.512821 312.307692 x3 = 20 8000 Constrói-se o polinómio: 92 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA p3(x) = 2310+(x−3)390+(x−3)(x−5)8.75+(x−3)(x−5)(x−7)×(−0.956637) = 5245.803167 Para a velocidade no instante t = 10 segundos tem-se f(10) ≈ p3(10) = 5245.803167. b) Para calcular o erro de truncatura, utiliza-se a expressão (A.13): |R3(x)| ≤ |(x − x0)(x − x1)(x − x2)(x − x3)||dd4| Para calcular a diferença dividida de ordem 4, é necessário incluir mais um ponto, o mais próximo possível - usa-se o ponto que tinha sido deixado de parte na alínea anterior. xi fi dd1 dd2 dd3 dd4 1 800 755 3 2310 −91.25 390 16.666667 5 3090 8.75 −0.927542 425 −0.956637 7 3940 −7.512821 312.307692 20 8000 Como tal, |R3(10)| ≤ |(10 − 3)(10 − 5)(10 − 7)(10 − 20)|| − 0.927542| = 973.9191 CONUM: DIFDIV Interpolação Polinomial - Diferenças Divididas Pontos Valores da Função 3.000000 2310.000000 5.000000 3090.000000 7.000000 3940.000000 20.000000 8000.000000 Ponto Interpolador = 10.000000 Grau do Polinómio = 3 Interpolação Directa 93 Tabela das Diferenças Diferenças de 1ł Ordem 390.000000 425.000000 312.307692 Diferenças de 2ł Ordem 8.750000 -7.512821 Diferenças de 3ł Ordem -0.956637 Aproximação ao Valor da Função no Ponto Interpolador: 5245.803167 94 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA 3. Considere a seguinte tabela da função f(x) que descreve a quantidade de pessoas que entram/saem num elevador em cada piso x numa subida entre o piso −2 e o piso 2: x −2 −1 0 1 2 f(x) a 2 1 0 4 a ∈ ℜ+. Determine o número de pessoas a que entram/saem no piso -2, de modo a que o polinómio interpolador de f(x) nos pontos da tabela dada seja de grau 3. Justifique. Resolução: Para o polinómo interpolador de Newton ser de grau 3, as diferenças divididas de ordem 3 devem ser iguais entre si e diferentes de zero. Calculando as diferenças divididas: xi fi dd1 dd2 dd3 −2 a 2 − a −1 2 −1 − (2 − a) 2 −1 −−1 − (2 − a) 2 3 0 1 0 −1 2.5 3 1 0 2.5 4 2 4 Para encontrar o valor de a, igualam-se as duas diferenças divididas de grau 3 (dd3), ou seja: −−1 − (2 − a) 2 3 = 2.5 3 ⇔ a = −2 Entram (saem) -2 pessoas no piso -2 (f(−2) = −2). 95 4. Considere a seguinte tabela que representa a deflexão em cm duma prancha de saltos, num salto de um atleta olímpico, em vários instantes de tempo de preparação: t 0 0.5 1.5 2 4 6 f(t) a 1.75 -1.25 -2 0 b a) Calcule os valores de a e b de forma a que a deflexão f(t) corresponda a um polinómio de grau ≤ 2. b) Construa o polinómio. Resolução: Mudança de variável: t → x. a) Para o polinómio ser de grau dois é necessário que todas as diferenças divididas de grau dois sejam iguais entre si e diferentes de zero. Calculam-se as diferenças divididas: xi fi dd1 dd2 0 a 1.75 − a 0.5 0.5 1.75 −3 − 1.75 − a 0.5 1.5 −3 1.5 −1.25 1 −1.5 2 −2 1 1 4 0 b 2 − 1 4 b 2 6 b As dd2 têm de ser todas iguais a 1, logo igualando as diferenças divididas desconhe- cidas a 1, encontram-se os valores de a e b. 96 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA −3 − 1.75 − a 0.5 1.5 = 1 ⇔ −3 − (1.75 − a 0.5 ) = 1.5 ⇔ 1.75 − a 0.5 = −4.5 ⇔ 1.75 − a = −2.25 ⇔ a = 4 b 2 − 1 4 ⇔ b 2 − 1 = 4 ⇔ b = 10 b) Para construir o polinómio recorre-se a (A.12). Como a tabela representa uma função f(x) que é um polinómio de grau igual ou inferior a dois, para a = 4 e b = 10, quaisquer 3 pontos podem ser usados para construir esse polinómio. Considere-se x0 = 0.5, x1 = 1.5 e x2 = 2. O polinómio é: p2(x) = 1.75 + (x − 0.5)(−3) + (x − 0.5)(x − 1.5)(1) 97 5. Um braço de um robô deve passar nos instantes t0, t1, t2, t3, t4 e t5 por posições pré- definidas θ(t0), θ(t1), θ(t2), θ(t3), θ(t4) e θ(t5), onde θ(t) é o ângulo (em radianos) que o braço do robô faz com o plano XOY. θ Z X Y ti 1 2 3 4 5 6 θi = θ(ti) 1 1.25 1.75 2.25 3 3.15 a) Com base nos dados da tabela, aproxime a trajectória do robô por uma spline cúbica completa. Indique também uma aproximação da posição do robô no instante t = 1.5. b) Calcule uma aproximação à velocidade do robô no instante t = 1.5. c) Calcule um limite superior do erro de truncatura que se comete quando se usa a derivada da spline calculada para aproximar a velocidade do robô. Resolução: Mudança de variável: t → x e θ → f. a) Não se conhece a expressão analítica de f(x) e como a spline é completa vai ser necessário estimar o valor das derivadas nos valores extremos. Assim, reservam-se dois pontos A e B, que não serão nós da spline, para estimar essas derivadas. Os nós da spline são: x0 = 1, x1 = 3, x2 = 4 e x3 = 6 (A = 2 e B = 5). f ′ 0 ≈ f(A) − f0 A − x0 = 1.25 − 1 2 − 1 = 0.25 f ′ n ≈ fn − f(B) xn − B = 3.15 − 3 6 − 5 = 0.15 Utilizando as expressões (A.17), determinam-se as equações respeitantes aos nós de fronteira, x0 e x3, respectivamente. • 2(x1 − x0)M0 + (x1 − x0)M1 = 6 x1 − x0 (f1 − f0) − 6f ′(x0) ⇔ 2(3 − 1)M0 + (3 − 1)M1 = 6 3 − 1(1.75 − 1) − 6 × 0.25 ⇔ 4M0 + 2M1 = 0.75 • 2(x3 − x2)M3 + (x3 − x2)M2 = 6f ′(x3) − 6 x3 − x2 (f3 − f2) ⇔ 2(6 − 4)M3 + (6 − 4)M2 = 6 × 0.15 − 6 6 − 4(3.15 − 2.25) ⇔ 4M3 + 2M2 = −1.8 Para os nós interiores, x1 (i = 1) e x2 (i = 2), respectivamente, utiliza-se (A.15): 102 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA O ponto x = 1.75 pertence ao 1° segmento ([1.5, 2.0]), logo deve utilizar-se i = 1 na expressão (A.14): s_3^1(x) = \frac{M_0}{6(x_1− x_0)}(x_1− x_0)^3 + \frac{M_1}{6(x_1− x_0)}(x− x_0)^3 + \left[\frac{f_0}{x_1− x_0}− \frac{M_0(x_1− x_0)}{6}\right](x− x_0) s_3^1(x) =\frac{7.3956565}{6(2.0− 1.5)} (1.75−1.5)^3 + \frac{4.9}{2.0− 1.5} (2.0−1.75)\ns_3^1(1.75) = 3.983424 \approx f(1.75) Um fio de diâmetro 1.75 tem uma resistência de aproximadamente 3.983424. CONUM: SPLINE Interpolação Polinomial - Splines Pontos Valores da Função 1.500000 4.900000 2.000000 3.300000 2.200000 3.000000 3.000000 2.000000 3.800000 1.750000 4.000000 1.500000 Ponto Interpolador = 1.750000 Spline Natural Polinômio do Segmento 1 da Spline 0.000000*pot(2.000000 − x, 3) + 2.486956*pot(x − 1.500000, 3) + 9.800000*(2.000000 − x) + 5.978261*(x − 1.500000) Valor da Spline no Ponto Interpolador = 3.983424 99 c) O erro da aproximação à derivada é dado por uma expressão em (A.19): |f ′(x)−s′ 3(x)| ≤ 1 24h3M4, em que M4 é o majorante da 4a derivada, i.e., é a diferença dividida de quarta ordem de maior módulo multiplicada por 4! (A.11) e h é a máxima distância entre 2 nós consecutivos da spline. Para uma melhor estimação do erro, incluem-se na tabela das diferenças divididas todos os pontos conhecidos, indiferentemente de serem ou não nós da spline: Calculando as diferenças divididas: xi fi dd1 dd2 dd3 dd4 1 1 0.25 2 1.25 0.125 0.5 −0.041667 3 1.75 0 0.020834 0.5 0.041667 4 2.25 0.125 −0.045834 0.75 −0.141667 5 3 −0.3 0.15 6 3.15 Erro: |f ′(x) − s′ 3(x)| ≤ 1 24 × 23 × | − 0.045834| × 4! = 0.366672 CONUM: SPLINE Interpolação Polinomial - Splines Pontos Valores da Função 1.000000 1.000000 3.000000 1.750000 4.000000 2.250000 6.000000 3.150000 Ponto Interpolador = 1.500000 100 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA Spline Completa Pontos Auxiliares xa = 2.000000 fa = 1.250000 xb = 5.000000 fb = 3.000000 Polinómio do Segmento 1 da Spline 0.013411*pot(3.000000 - x, 3) + 0.004427*pot(x - 1.000000, 3) + 0.446354*(3.000000 - x) + 0.857292*(x - 1.000000) Valor da Spline no Ponto Interpolador = 1.143994 MATLAB: >> x=[1 2 3 4 5 6] >> y=[1 1.25 1.75 2.25 3 3.15] >> sol=spline(x,y,1.5) sol = 1.0652 >> res=spline(x,y) >> res.coefs ans = -0.0761 0.3533 -0.0272 1.0000 -0.0761 0.1250 0.4511 1.2500 0.1306 -0.1033 0.4728 1.7500 -0.1961 0.2883 0.6578 2.2500 -0.1961 -0.3000 0.6461 3.0000 Expressão do polinómio no 1o segmento [1, 2]: s1 3(x) = −0.0761(x − 1)3 + 0.3533(x − 1)2 − 0.0272(x − 1) + 1.0000. Expressão do polinómio no 2o segmento [2, 3]: s2 3(x) = −0.0761(x − 2)3 + 0.1250(x − 2)2 + 0.4511(x − 2) + 1.2500. ... 101 6. A resistência de um certo fio de metal, f(x), varia com o diâmetro desse fio, x. Foram medidas as resistências de 5 fios de diversos diâmetros: xi 1.5 2.0 2.2 3.0 3.8 f(xi) 4.9 3.3 3.0 2.0 1.75 Use uma spline cúbica natural para estimar a resistência de um fio de diâmetro 1.75. Resolução: Existem 5 nós (x0 = 1.5, x1 = 2.0, x2 = 2.2, x3 = 3.0 e x4 = 3.8) e por isso 4 segmentos. Inicia-se a resolução pelo cálculo das equações para os 3 nós interiores (A.15): i = 1 (x1 − x0)M0 + 2(x2 − x0)M1 + (x2 − x1)M2 = 6 x2 − x1 (f2 − f1) − 6 x1 − x0 (f1 − f0) ⇔ (2.0 − 1.5)M0 + 2(2.2 − 1.5)M1 + (2.2 − 2.0)M2 = 6 2.2 − 2.0(3.0 − 3.3) − 6 2.0 − 1.5(3.3 − 4.9) ⇔ 0.5M0 + 1.4M1 + 0.2M2 = 10.2 i = 2 (x2 − x1)M1 + 2(x3 − x1)M2 + (x3 − x2)M3 = 6 x3 − x2 (f3 − f2) − 6 x2 − x1 (f2 − f1) ⇔ (2.2 − 2.0)M1 + 2(3.0 − 2.0)M2 + (3.0 − 2.2)M3 = 6 3.0 − 2.2(2.0 − 3.0) − 6 2.2 − 2.0(3.0 − 3.3) ⇔ 0.2M1 + 2M2 + 0.8M3 = 1.5 i = 3 (x3 − x2)M2 + 2(x4 − x2)M3 + (x4 − x3)M4 = 6 x4 − x3 (f4 − f3) − 6 x3 − x2 (f3 − f2) ⇔ (3.0 − 2.2)M2 + 2(3.8 − 2.2)M3 + (3.8 − 3.0)M4 = 6 3.8 − 3.0(1.75 − 2.0) − 6 3.0 − 2.2(2.0 − 3.0) ⇔ 0.8M2 + 3.2M3 + 0.8M4 = 5.625 O resultado é um sistema de três equações e cinco incógnitas. No entanto, como se trata de uma spline cúbica natural, a concavidade é nula nos nós de fronteira, i.e., M0 = 0 e M4 = 0, originando um sistema com três equações e três incógnitas:      1.4 0.2 0 | 10.2 0.2 2 0.8 | 1.5 0 0.8 3.2 | 5.625      Resolvendo o sistema linear por EGPP, obtêm-se os valores de M1 = 7.395665,M2 = −0.769657 e M3 = 1.950227. 103 7. A função s3(t) representa uma spline cúbica que descreve a quantidade de um nu- triente na corrente sanguínea na 1a hora (0 ≤ t ≤ 1) após a administração dum complexo vitamínico e entre a 2a e 5a horas (1 ≤ t ≤ 5): s3(t) =    s1 3(t) = at3 + 2.9t + 1, 0 ≤ t ≤ 1 s2 3(t) = −0.025t3 + bt2 + 2.525t + 1.125, 1 ≤ t ≤ 5 Calcule os valores reais de a e b para os quais a spline é natural para aproximar uma função nos pontos t = 0, 1 e 5. Resolução: Mudança de variável: t → x. Para a spline ser natural, os valores das segundas derivadas da função nos extremos devem ser iguais a zero. Ao mesmo tempo, deve assegurar-se que os valores da função em t = 1, único ponto interior, devem ser idênticos, bem como os valores da primeira e segunda derivadas (continuidade até à 2a ordem). Assim, determina-se a expressão da primeira derivada: s′ 3(x) =    s1′ 3 (x) = 3ax2 + 2.9, 0 ≤ x ≤ 1 s2′ 3 (x) = −0.075x2 + 2bx + 2.525, 1 ≤ x ≤ 5 e da segunda derivada: s′′ 3(x) =    s1′′ 3 (x) = 6ax, 0 ≤ x ≤ 1 s2′′ 3 (x) = −0.15x + 2b, 1 ≤ x ≤ 5 Sabendo que a spline é natural, s1′′ 3 (0) = 0 e s2′′ 3 (5) = 0. s1′′ 3 (0) = 6a × 0 = 0 ⇔ 0a = 0, ou seja, o valor de a pode ser qualquer valor real. s2′′ 3 (5) = −0.15 × 5 + 2b = 0 ⇔ b = 0.375. O valor de a pode ser qualquer valor real. No entanto, sabe-se que, pela continuidade até à 2a ordem, s1 3(1) = s2 3(1), s1′ 3 (1) = s2′ 3 (1) e s1′′ 3 (1) = s2′′ 3 (1). Assim, substituindo b 104 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA pelo valor encontrado, b = 0.375, nas equações anteriores, tem-se: s1 3(1) = s2 3(1) = 4, s1′ 3 (1) = s2′ 3 (1) = 3.2 e s1′′ 3 (1) = s2′′ 3 (1) = 0.6 Substituindo na expressão de s2′′ 3 (1): 6ax = 0.6 ⇔ 6a × 1 = 0.6 ⇔ a = 0.1 Assim, a = 0.1 e b = 0.375 105 8. A distância requerida para parar um automobilista é função da velocidade a que ele se desloca. Os seguintes dados experimentais foram recolhidos para quantificar essa relação: vel (Km/h) 15 20 25 30 40 50 distância (m) 16 20 34 40 60 90 Estime a distância necessária para parar um carro que se desloca a uma velocidade de 45 Km/h, utilizando uma spline cúbica completa Resolução: Mudança de variável: vel → x e distância → f(x). Não se conhece a expressão analítica de f(x) e como a spline é completa vai ser necessário estimar o valor das derivadas nos valores extremos. Assim, reservam-se dois pontos A e B (sendo A (x = 20) o ponto mais próximo e superior a x0, e B (x = 40) o ponto mais próximo e inferior a x3), que não serão nós da spline, para estimar essas derivadas. Os nós da spline são: x0 = 15, x1 = 25, x2 = 30 e x3 = 50. Assim, tem-se: f ′ 0 ≈ f(A) − f0 A − x0 = 20 − 16 20 − 15 = 0.8 e f ′ n ≈ fn − f(B) xn − B = 90 − 60 50 − 40 = 3. Os pontos A e B são excluídos dos nós da spline e utilizam-se apenas os valores: vel (Km/h) 15 25 30 50 distância (m) 16 34 40 90 Determinam-se as equações respeitantes aos nós da fronteira (x0 e x3, respectiva- mente), utilizando a expressão (A.17). • 2(x1 − x0)M0 + (x1 − x0)M1 = 6 x1 − x0 (f1 − f0) − 6f ′(x0) ⇔ 2(25 − 15)M0 + (25 − 15)M1 = 6 25 − 15(34 − 16) − 6 × 0.8 ⇔ 20M0 + 10M1 = 6 106 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA • 2(x3 − x2)M3 + (x3 − x2)M2 = 6f (x3) − 6 x3 − x2 (f3 − f2) ⇔ 2(50 − 30)M3 + (50 − 30)M2 = 6 × 3 − 6 50 − 30 (90 − 40) ⇔ 40M3 + 20M2 = 3 Para os nós interiores, utiliza-se (A.15), com i = 1 e i = 2, respectivamente: • (x1−x0)M0 +2(x2−x0)M1+(x2−x1)M2 = 6 x2 − x1 (f2−f1)− 6 x1 − x0(f1−f0) ⇔ (25−15)M0+ 2(30−15)M1 +(30−25)M2 = 6 30 − 25 (40− 34)− 6 25 − 15 (34− 16) ⇔ 10M0 + 30M1 + 5M2 = −3.6 • (x2−x1)M1 +2(x3−x1)M2+(x3−x2)M3 = 6 x3 − x2 (f3−f2)− 6 x2 − x1 (f2−f1) ⇔ (30−25)M1+ 2(50−25)M2+(50−30)M3 = 6 50 − 30 (90−40)− 6 30 − 25 (40−34) ⇔ 5M1 + 50M2+20M3 = 7.8 Deste modo, é possível construir um sistema linear de 4 equações a 4 incógnitas: ⎡ ⎢⎢⎢⎢⎢⎢⎣ 20 10 0 0 0 0 0 5 0 0 0 0 5 20 10 0 0 0 20 40 ⎤ ⎥⎥⎥⎥⎥⎥⎦ ⎡ ⎢⎢⎢⎢⎢⎢⎣ 6 −3.6 7.8 3 ⎤ ⎥⎥⎥⎥⎥⎥⎦ Resolvendo o sistema linear por EGPP, obtém-se os valores M0 = 0.451538, M1 = −0.303077, M2 = 0.195385 e M3 = −0.022692. O ponto x = 45 para o qual se pretende calcular a função encontra-se situado no 3º segmento [30,50]. Utilizando a expressão (A.14) com i = 3, tem-se: s3 3(x) = M2 6(x3 − x2) (x3 − x)3 + M3 6(x3 − x2) (x − x2)3 + f2 x3 − x2 − M2(x3 − x2) 6 (x3 − x)+ + ⎡ ⎢⎢⎢⎣ f3 x3 − x2 − M3(x3 − x2) 6 ⎤ ⎥⎥⎥⎦ (x − x2) s3 3(x) = 0.195385 6(50−30) (50− x)3 + −0.022692 6(50 − 30) (x − 30)3 + ⎡ ⎢⎢⎢⎣ 40 50 − 30 − 0.195385(50 − 30) 6 ⎤ ⎥⎥⎥⎦ (50 − x) + ⎡ ⎢⎢⎢⎣ 90 50 − 30 − −0.022692(50 − 30) 6 ⎤ ⎥⎥⎥⎦ (x − 30) Para x = 45: s3 3(45) = 0.195385 6(50−30) (50−45)3+ −0.022692 6(50−30) (45−30)3+ ⎡ ⎢⎢⎢⎣ 40 50 − 30 − 0.195385(50 − 30) 6 ⎤ ⎥⎥⎥⎦ (50 − 45)+ ⎡ ⎢⎢⎢⎣ 90 50 − 30 − −0.022692(50 − 30) 6 ⎤ ⎥⎥⎥⎦ (45 − 30) = 74.943510 ≈ f(45) 107 CONUM: SPLINE Interpolação Polinomial - Splines Pontos Valores da Função 15.000000 16.000000 25.000000 34.000000 30.000000 40.000000 50.000000 90.000000 Ponto Interpolador = 45.000000 Spline Completa Pontos Auxiliares xa = 20.000000 fa = 20.000000 xb = 40.000000 fb = 60.000000 Polinómio do Segmento 3 da Spline 0.001628*pot(50.000000 - x, 3) + -0.000189*pot(x - 30.000000, 3) + 1.348718*(50.000000 - x) + 4.575641*(x - 30.000000) Valor da Spline no Ponto Interpolador = 74.943510 108 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA 9. Considere a função definida por x −2 0 1 2 f(x) −8 0 1 8 que corresponde às coordenadas (x, f(x)) de pontos num ecrã de computador para processamento gráfico. Sabendo que s1′′ 3 (−2) = 12 e sn′′ 3 (2) = 20 estime o valor de f(−1) através de uma spline cúbica. Resolução: O valor de s1′′ 3 (−2) e sn′′ 3 (2) corresponde a M0 = 12 e Mn = 20, respectivamente. O facto de serem valores diferentes de zero implica que a spline é completa. Os nós serão x0 = −2, x1 = 0, x2 = 1 e x3 = 2. São 2 os pontos interiores (A.15): i = 1 (x1 − x0)M0 + 2(x2 − x0)M1 + (x2 − x1)M2 = 6 x2 − x1 (f2 − f1) − 6 x1 − x0 (f1 − f0) ⇔ (0 − (−2))M0 + 2(1 − (−2))M1 + (1 − 0)M2 = 6 1 − 0(1 − 0) − 6 0 − (−2)(0 − (−8)) ⇔ 2M0 + 6M1 + 1M2 = −18 i = 2 (x2 − x1)M1 + 2(x3 − x1)M2 + (x3 − x2)M3 = 6 x3 − x2 (f3 − f2) − 6 x2 − x1 (f2 − f1) ⇔ (1 − 0)M1 + 2(2 − 0)M2 + (2 − 1)M3 = 6 2 − 1(8 − 1) − 6 1 − 0(1 − 0) ⇔ 1M1 + 4M2 + 1M3 = 36 As duas equações resultantes, embora com 4 incógnitas (M0, M1, M2 e M3), po- dem ser simplificadas, uma vez que os valores de M0 e M3 são dados no enunciado, originando o sistema linear seguinte.    2M0 + 6M1 + 1M2 = −18 1M1 + 4M2 + 1M3 = 36 ⇔    2 × 12 + 6M1 + 1M2 = −18 1M1 + 4M2 + 1 × 20 = 36 ⇔    6M1 + 1M2 = −42 1M1 + 4M2 = 16 Resolvendo o sistema linear por EGPP, obtêm-se os valores M1 = −8 e M2 = 6. O ponto interpolador é x = −1, que se encontra no primeiro segmento [-2,0]. Utilizando a expressão (A.14) com i = 1, obtém-se: 109 s1 3(x) = M0 6(x1 − x0) (x1 − x)3 + M1 6(x1 − x0) (x − x0)3 + ⎡ ⎢⎢⎢⎣ f0 x1 − x0 − M0(x1 − x0) 6 ⎤ ⎥⎥⎥⎦ (x1 − x)+ + ⎡ ⎢⎢⎢⎣ f1 x1 − x0 − M1(x1 − x0) 6 ⎤ ⎥⎥⎥⎦ (x − x0) s1 3(x) = 12 6(0−(−2)) (−2 − x)3 + −8 6(0−(−2)) (x − (−2))3 + ⎡ ⎢⎢⎢⎣ −8 0 − (−2) − −8(0−(−2)) 6 ⎤ ⎥⎥⎥⎦ (x − (−2))+ ⎡ ⎢⎢⎢⎣ 12 0−(−2) (0−(−2)) 6 ⎤ ⎥⎥⎥⎦ (0 − x)+ Para x = −1: s1 3(−1) = −5 ≈ f(−1) 110 CAPÍTULO 5. INTERPOLAÇÃO NUMÉRICA 10. A seguinte função segmentada s3(x) no intervalo [0, 3], representa o lucro obtido na venda de um produto sazonal. No 1o mês de vendas o lucro é representado por s1 3(x) e nos 2o e 3o meses é descrito por s2 3(x). Poderá a função segmentada s3(x) representar uma spline cúbica? Justifique. s3(x) =    s1 3(x) = 3 x3 − x2 + x − 2, 0 ≤ x ≤ 1 s2 3(x) = 2 x3 + 2x − 3, 1 ≤ x ≤ 3 Resolução: Para esta função cúbica segmentada representar uma spline, deve ser continuamente diferenciável até à segunda ordem. Nos pontos interiores (x = 1) deverá verificar-se: s1 3(1) = s2 3(1), s1′ 3 (1) = s2′ 3 (1), s1′′ 3 (1) = s2′′ 3 (1). A função tem o mesmo valor para x = 1 nos dois segmentos: s1 3(1) = 1 e s2 3(1) = 1. Analisando a primeira derivada: s′ 3(x) =    s1′ 3 (x) = 9x2 − 2x + 1, 0 ≤ x ≤ 1 s2′ 3 (x) = 6x2 + 2, 1 ≤ x ≤ 3 s1′ 3 (1) = 8 e s2′ 3 (1) = 8. A primeira derivada da função tem o mesmo valor nos dois segmentos, para x = 1. Verificando a segunda derivada: s′′ 3(x) =    s1′′ 3 (x) = 18x − 2, 0 ≤ x ≤ 1 s2′′ 3 (x) = 12x, 1 ≤ x ≤ 3 s1′′ 3 (1) = 16 e s2′′ 3 (1) = 12. A segunda derivada da função tem valores diferentes para x = 1 em cada um dos segmentos. Este facto indica que a função s3(x) não pode representar uma spline. Capítulo 6 Mínimos Quadrados Dada uma função definida num intervalo [a, b] ∈ ℜ por uma tabela matemática com m pontos xi x1 x2 . . . xm fi f1 f2 . . . fm em que a ≤ x1 < x2 < . . . < xm ≤ b, ou por uma relação funcional f(x), pretende-se calcular um modelo linear nos parâmetros, M(x), que reflicta na generalidade, o compor- tamento dos dados. Esse modelo é construído com o objectivo de minimizar S (A.44), o somatório do quadrado dos resíduos: m S = ∑ [f(xi) − M(xi)]2 . i=1 O modelo M(x) pode ser um polinómio pn(x) de grau menor ou igual a n (A.38). Esse polinómio vai ser construído através duma relação de recorrência (A.39), (A.40) e (A.41) que gera polinómios ortogonais. O modelo pode também ser não polinomial (A.42) obri- gando à resolução dum sistema de equações normais (A.43) para o cálculo dos parâmetros. O resíduo mede a proximidade de f(x) em relação a M(x), ie, funciona como um mecanismo de avaliação do modelo M(x). Vão ser calculados modelos lineares nos parâmetros, do tipo polinomial e não polino- mial. 111 112 CAPÍTULO 6. MÍNIMOS QUADRADOS 1. Na realização de uma experiência, determinaram-se os seguintes valores da capacidade de calor C a várias temperaturas T, de um metal: \begin{tabular}{|c|c|c|c|c|c|c|} \hline T & -50 & -20 & 10 & 70 & 100 & 120 \\ \hline C & 0.125 & 0.128 & 0.134 & 0.144 & 0.150 & 0.155 \\ \hline \end{tabular} Use a técnica dos mínimos quadrados para calcular: a) um modelo linear que estime C como função de T. b) um modelo quadrático que faça prever a relação entre C e T. c) qual dos modelos lhe parece mais adequado? Justifique. Resolução: Mudança de variável: T → x e C → f(x). a) O cálculo de um modelo linear envolve a criação de um polinômio de primeiro grau, com o auxílio da expressão p1(x) = c0P0(x) + c1P1(x) (A.38). Primeiro calculam-se os P’s usando a relação de recorrência (A.39), em que P−1 = 0 e P0 = 1. Utiliza-se também (A.40) para encontrar os valores de Bi e Ci necessários no cálculo de Pi. A construção de uma tabela com os somatórios simplifica o exercício. \begin{tabular}{|c|c|c|c|c|c|} \hline x & f(x) & P0 & P1 & P12 & fP1 \\ \hline -50 & 0.125 & 1 & -88.333333 & 7802.777779 & -11.041667 \\ -20 & 0.128 & 1 & -58.333333 & 3402.777739 & -7.466667 \\ 10 & 0.134 & 1 & -21.333333 & 802.777759 & -3.796667 \\ 70 & 0.144 & 1 & 31.666667 & 1002.777799 & 4.560000 \\ 100 & 0.150 & 1 & 61.666667 & 3802.777819 & 9.250000 \\ 120 & 0.155 & 1 & 81.666667 & 6669.444499 & 12.683333 \\ \hline ∑ 230 & 0.836 & 6 & - & 23483.333333 & 4.163334 113 CAPÍTULO 6. MÍNIMOS QUADRADOS B0 = ∑xP20∑P20 = 2306 = 38.333333 P1(x) = (x−B0)P0(x)−C0P−1(x) ⇒ P1(x) = x−x−B0 = x−38.333333 Após o cálculo dos P’s, é necessário calcular os c’s (A.41): c0 = ∑fP0∑P20 = 0.8366 = 0.139333 c1 = ∑fP1∑P21 = 4.16333423483.333333 = 0.000177 O polinômio p1(x) tem então a seguinte expressão: p1(x) = c0P0(x) + c1P1(x) ⇒ p1(x) = 0.139333 + 0.000177(x−38.333333) b) Para calcular um modelo quadrático recorre-se a um polinômio do tipo: p2(x) = c0P0(x) + c1P1(x) + c2P2(x) Como já tinham sido calculados na alínea a) os termos c0P0 e c1P1, apenas se tem de calcular P2(x) (i = 1 em (A.39)). P2(x) = (x−B1)P1(x)−C1P0(x) ⇒ P2(x) = (x−B1)(x−38.333333)−C1 \begin{tabular}{|c|c|c|c|c|c|c|} \hline x & f(x) & xP21 & P2 & P22 & fP2 \\ \hline -50 & 0.125 & -390138.885944 & 3514.407438 & 12351059.640050 & 439.300930 \\ -20 & 0.128 & -68055.554778 & -758.410182 & 575186.004209 & -97.076503 \\ 10 & 0.134 & 8027.777589 & -3231.227802 & 10440833.108620 & -432.984525 \\ 70 & 0.144 & 70194.445922 & -2776.863042 & 7710968.354199 & -399.868278 \\ 100 & 0.150 & 380277.781889 & 150.319338 & 22595.903637 & 22.547901 \\ 120 & 0.155 & 800333.339867 & 3101.774258 & 9621003.547397 & 480.775010 \\ \hline ∑ - & - & 800638.904544 & - & 40721646.557843 & 12.694534 \\ \hline \end{tabular} B1 = ∑xP21∑P21 = 800638.90454423483.333392 = 34.093921 114 CAPÍTULO 6. MÍNIMOS QUADRADOS C1 = ∑P21∑P20 = 23483.3333926 = 3913.888889 P2(x) = (x−34.093921)(x−38.333333)−3913.888889 c2 = ∑fP2∑P22 = 12.69453440721646.557843 = 0.000000 O modelo quadrático (polinômio de segundo grau) tem a seguinte expressão: p2(x) = 0.139333 × 1 + 0.000177(x−38.333333) + 0.000000 × [(x−34.093921)(x−38.333333)−3913.888889] ⇔ p2(x) = 0.139333 + 0.000177(x−38.333333) c) Neste caso, os modelos linear e quadrático são idênticos, uma vez que o coeficiente c2 é nulo e p2(x) reduz-se ao polinômio p1(x). CONUM: MQPOLI Mínimos Quadrados - Modelo Polinomial Linear Pontos Valores da Função -50.000000 0.125000 -20.000000 0.128000 10.000000 0.134000 70.000000 0.144000 100.000000 0.150000 120.000000 0.155000 Ponto Interpolador = 40.000000 Coeficientes B e C da Relação de Recorrência B[0] = 38.333333 Coeficientes do Polinômio c[0] = 0.139333 c[1] = 0.000177 Soma dos Quadrados dos Resíduos = 0.000005 Valor do Polinômio no Ponto Interpolador = 0.139629 ================================================================================ MQPOLI Mínimos Quadrados - Modelo Polinomial Linear Pontos Valores da Função -50.000000 0.125000 115 -20.000000 0.128000 10.000000 0.134000 70.000000 0.144000 100.000000 0.150000 120.000000 0.155000 Ponto Interpolador = 40.000000 Coeficientes B e C da Relação de Recorrência B[0] = 38.333333 B[1] = 34.093920 C[1] = 3913.888889 Coeficientes do Polinómio c[0] = 0.139333 c[1] = 0.000177 c[2] = 0.000000 Soma dos Quadrados dos Resíduos = 0.000001 Valor do Polinómio no Ponto Interpolador = 0.138412 MATLAB: >> x=[-50 -20 10 70 100 120] >> y=[0.125 0.128 0.134 0.144 0.150 0.155] >> [P,s]=polyfit(x,y,1) P = 0.000177 0.132537 s = R: [2x2 double] df: 4 normr: 0.002285 p1(x) = 0.000177x + 0.132537. O somatório do quadrado do erro é 0.0022852 = 0.000005. MATLAB: >> [P,s]=polyfit(x,y,2) P = 0.000000 0.000155 0.131725 s = R: [3x3 double] df: 3 normr: 0.001124 p2(x) = 0.000000x2 + 0.000155x + 0.131725. O somatório do quadrado do erro é 0.0011242 = 0.000001. 116 CAPÍTULO 6. MÍNIMOS QUADRADOS 2. O preço de um veículo usado é função de diversos factores, entre os quais figuram os seguintes: o modelo do veículo em causa, o preço do veículo novo, a idade do veículo usado, o seu estado de conservação e a relação procura/oferta no mercado de veículos em segunda mão. Considerando apenas o caso de veículos devidamente conservados e supondo que o factor procura/oferta do mercado em segunda mão não variou significativamente, o preço de um veículo usado de determinado modelo pode ser explicado através de uma relação do tipo f(I) = PV U PV N em que PV U é o preço do veículo usado, PV N é o preço do veículo novo e I é a idade do veículo usado. Na tabela seguinte apresentam-se valores de f e I para veículos de um determinado modelo, que se supõem obtidos nas condições acima mencionadas. I (anos) 1 2 3 4 5 6 f(I) 0.843 0.753 0.580 0.520 0.452 0.414 Determine o modelo do tipo M(I; c1, c2) = c1+c2e−0.1I que melhor se ajusta à função f(I) no sentido dos mínimos quadrados. Nestes cálculos utilize apenas a informação relativa aos 3 pontos (da tabela) que correspondem a I = 1, 3, 6. Resolução: Mudança de variável: I → x e f(I) → f(x). Para encontrar o modelo pretendido que é do tipo não polinomial, M(x) = c1φ1(x)+ c2φ2(x), devem calcular-se c1 e c2. O modelo tem a forma apresentada na expressão (A.42), com φ1 = 1 e φ2 = e−0.1x. Para determinar os coeficientes do modelo, é necessária a resolução do sistema linear indicado em (A.43), cujos valores são obtidos através dos cálculos apresentados na seguinte tabela: x f(x) φ₁ φ₂ φ₁² φ₂² φ₁φ₂ fφ₁ fφ₂ 1 0.843 1 0.904837 1 0.818731 0.904837 0.843 0.762778 3 0.580 1 0.740818 1 0.548812 0.740818 0.580 0.429675 6 0.414 1 0.548812 1 0.301194 0.548812 0.414 0.227208 Σ - - - 3 1.668737 2.194467 1.837 1.419661 Constroi-se o sistema linear em c₁ e c₂: ⎛ 3 2.194467 ⎞ ⎛ c₁ ⎞ = ⎛ 1.837 ⎞ ⎝ 2.194467 1.668737 ⎠ ⎝ c₂ ⎠ ⎝ 1.419661 ⎠ Resolvendo o sistema linear por EGPP obtêm-se os valores dos parâmetros c₁ = -0.262061 e c₂ = 1.195362. O modelo obtido é M(x) = -0.262061 + 1.195362e⁻⁰.¹ˣ CONUM: MQLINE Mínimos Quadrados - Modelo Não Polinomial Linear Vector de Funções 1 exp(-0.1*x1) Pontos Valores da Função 1.000000 0.843000 3.000000 0.580000 6.000000 0.414000 Ponto Interpolador = 2.000000 Multiplicadores da Eliminação de Gauss m[2,1] = -0.731489 Sistema possível e determinado. Solução única. Parâmetros do Modelo c[1] = -0.262065 c[2] = 1.195367 Soma dos Quadrados dos Resíduos = 0.002842 Valor do Modelo no Ponto Interpolador = 0.716619 CAPÍTULO 6. MÍNIMOS QUADRADOS 3. O consumo de gás natural sofre uma redução significativa durante os meses de Verão. Na tabela seguinte estão registados alguns valores recolhidos ao longo do ano de 2006. mes 1 3 4 6 9 12 consumo 20.0 7.5 6.5 7.0 10.0 A Uma companhia de gás sugeriu um modelo do tipo M(x; c₁, c₂) = c₁x² + c₂ 1/x para estimar o consumo de gás em qualquer altura do ano. No sentido dos mínimos quadrados e considerando a amostra de 6 pontos, (a) comecar por apresentar o sistema de equações lineares que deve construir para calcular os parâmetros c₁ e c₂, em função de A; (b) considerando A = 15.0 apresente o modelo sugerido. Mudança de variável: mês → x e consumo de gás → f(x). a) O modelo não polinomial disponibilizado é do tipo da equação (A.42), com φ₁ = x² e φ₂ = 1/x. Para obter os coeficientes do modelo são necessários cálculos apresentados na seguinte tabela. x f(x) φ₁ φ₂ φ₁² φ₂² φ₁φ₂ fφ₁ fφ₂ 1 20.0 1 1 1 1 1 20 20 3 7.5 9 0.333333 81 0.111111 3 67.5 2.5 4 6.5 16 0.25 256 0.0625 4 104 1.625 6 7.0 36 0.166667 1296 0.027778 6 252 1.166667 9 10.0 81 0.111111 6561 0.012346 9 810 1.111111 12 A 144 0.083333 20736 0.006944 12 144A 0.083333A ∑ - - - 28931 1.220679 35 1253.5+144A 26.402778+0.083333A 119 O sistema linear para a determinação dos coeficientes é dado pela expressão (A.43):   28931 35 35 1.220679     c1 c2   =   1235.5 + 144A 26.402778 + 0.083333A   b) Com A = 15, o sistema toma a forma:   28931 35 35 1.220679     c1 c2   =   3413.5 27.652778   A resolução do sistema linear por EGPP dá c1 = 0.093837 e c2 = 19.963061 O modelo é M(x) = 0.093837x2 + 19.963061 x CONUM: MQLINE Mínimos Quadrados - Modelo Não Polinomial Linear Vector de Funções pot(x1,2) 1/x1 Pontos Valores da Função 1.000000 20.000000 3.000000 7.500000 4.000000 6.500000 6.000000 7.000000 9.000000 10.000000 12.000000 15.000000 Ponto Interpolador = 5.000000 Parâmetros do Modelo c[1] = 0.093837 c[2] = 19.963061 Soma dos Quadrados dos Resíduos = 0.153952 Valor do Modelo no Ponto Interpolador = 6.338533 121 A radiação média no mês de Agosto é dada pelo valor da função com x = 8: M(8) = 24.720339 × 8 + 148.314761sen(8) = 344.499144 A radiação média durante o mês de Agosto é aproximadamente 344.499144. CONUM: MQLINE Mínimos Quadrados - Modelo Não Polinomial Linear Vector de Funções x1 sen(x1) Pontos Valores da Função 1.000000 122.000000 3.000000 188.000000 6.000000 270.000000 10.000000 160.000000 12.000000 120.000000 Ponto Interpolador = 8.000000 Parâmetros do Modelo c[1] = 24.720341 c[2] = 148.314721 Soma dos Quadrados dos Resíduos = 45461.075610 Valor do Modelo no Ponto Interpolador = 344.499120 122 CAPÍTULO 6. MÍNIMOS QUADRADOS 5. O custo de investimento (C) em construção civil de um arejador num sistema de lamas activadas numa Estação de Tratamento de Águas Residuais depende do volume (v) do tanque da seguinte forma: C(v, c1, c2) = c1vc2 em que c1 e c2 são parâmetros a estimar pela técnica dos mínimos quadrados a partir dos dados recolhidos de uma construtora. vi em mil m3 0.4 0.6 1 1.3 Ci em milhares de euros 87 160 190 366 Estime os parâmetros c1 e c2 do modelo dado anteriormente, recorrendo à seguinte transformação que transforma o modelo dado num modelo polinomial de grau um: ln(C(v; c1; c2)) = ln(c1) + c2 ln(v) C = c1 + c2v Comece por calcular os parâmetros c1 e c2 do modelo polinomial usando a técnica dos mínimos quadrados, com base nos valores da tabela: vi = ln(vi) -0.916 -0.511 0 0.262 Ci = ln(Ci) 4.466 5.075 5.247 5.903 e posteriormente apresente os valores solicitados. Resolução: O modelo foi linearizado, sendo agora do tipo polinomial de grau 1 (A.38): p1(x) = c0P0(x) + c1P1(x) 132 CAPÍTULO 7. INTEGRAÇÃO NUMÉRICA b) Resolução Manual: Para calcular uma estimativa do erro de truncatura efectuado no cálculo, é necessário calcular os erros de truncatura provenientes de cada uma das fórmulas compostas de integração. Por isso: |eTotal| = |eCS(5)| + |eCT(2)| + |eCS(2.5)| |eCS(5)| = | − h4 180(b − a)f (iv)(η)| = | − 54 180(10 − 0) × 0| = 0 |eCT(2)| = | − h2 12(b − a)f ′′(η)| = | − 22 12(20 − 10) × 0| = 0 |eCS(5)| = | − h4 180(b − a)f (iv)(η)| = | − 2.54 180 (30 − 20) × 0| = 0 Ou seja, o erro de truncatura total é nulo, i.e., a estimativa do integral de 0 a 30 é exacta. 137 b) Para calcular uma estimativa do erro de truncatura no intervalo [5,15], recorre-se à expressão do erro (A.28): eCT = −h2 12(b − a)f ′′(η) No caso deste cálculo, apenas existem dois pontos, e são necessários pelo menos três para estimar a segunda derivada através das diferenças divididas de 2a ordem. Assim sendo, utiliza-se a ponto mais próximo do intervalo utilizado, que no caso corresponde a x = 2.5. Calculam-se as diferenças divididas: x F(x)cos(θ(x)) dd1 dd2 2.5 4.351270 −1.128626 5 1.529704 0.147673 0.717284 15 8.702540 Sabendo que f ′′ = dd2 × 2!, o valor do erro será: eCT ≤ | − 102 12 (15 − 5) × 0.147673 × 2!| = 24.612167 147 QUADRA Integração - Fórmulas de Newton-Cotes Pontos Valores da Função 0.000000 0.000000 0.500000 0.350000 1.000000 0.550000 1.500000 0.900000 Usou a Fórmula dos 3 Oitavos (h = 0.500000) Aproximação ao Integral = 0.675000 b) O erro de truncatura da fórmula composta do Trapézio é dado por (A.28), e deve ser inferior a 0.3, logo: |h2 12(b − a)f ′′(η)| ≤ 0.3 Deriva-se a expressão de f(x) para obter a segunda derivada da expressão. f(x) = 0.5x2 − 0.15x ; f ′(x) = x − 0.15 ; f ′′(x) = 1 A segunda derivada é constante e igual a 1, logo podem substituir-se os valores na expressão do erro: |h2 12(7.5 − 1.5) × 1| ≤ 0.3 ⇔ h ≤ 0.7746 Como h = b−a n , vem: 6 n < 0.7746 ⇔ n > 7.745933 Logo n = 8 e h = 7.5 − 1.5 8 = 0.75. xi 1.5 2.25 3 3.75 4.5 5.25 6 6.75 7.5 f(xi) 0.9 2.19375 4.05 6.46875 9.45 12.99375 17.1 21.76875 27 Utilizando a fórmula composta do Trapézio (A.27), obtém-se: v(x) ≈ 0.75 2 (0.9 + 2 × 2.19375+2×4.05+2×6.46875+2×9.45+2×12.99375+2×17.1+2×21.76875+27) = 65.98125 A velocidade do dragster é aproximadamente igual a 65.98125 m/s. CONUM: 148 CAPÍTULO 7. INTEGRAÇÃO NUMÉRICA QUADRA Integração - Fórmulas de Newton-Cotes Pontos Valores da Função 1.500000 0.900000 2.250000 2.193750 3.000000 4.050000 3.750000 6.468750 4.500000 9.450000 5.250000 12.993750 6.000000 17.100000 6.750000 21.768750 7.500000 27.000000 Usou a Fórmula de Simpson (h = 0.750000), Aproximação ao Integral = 65.700000 c) Na alínea a) foi utilizada a fórmula composta dos Três Oitavos (A.32). ET3/8 = | − h4 80(b − a)f (iv)(η)|, η ∈ [a, b] Para a estimação de f (iv) recorre-se ao cálculo de diferenças divididas de quarta ordem. São necessários 5 pontos, e apenas se dispõe de 4, pelo que se utiliza o ponto mais próximo do intervalo seleccionado. ti ai dd1 dd2 dd3 dd4 0 0 0.7 0.5 0.35 −0.3 0.4 0.4 1 0.55 0.3 −0.045714 0.7 0.297143 1.5 0.9 0.82 1.725 2.25 2.19375 Sabendo que f (iv) pode ser aproximado por dd4 × 4!, obtém-se: ET3/8 ≤ |h4 80 × (b − a)M4| = |0.54 80 (1.5 − 0) × 0.045714 × 4!| = 0.001286 O erro da fórmula composta dos Três Oitavos da alínea a) é de aproximadamente 0.001286. Capítulo 8 Diferenciação Numérica A dificuldade em encontrar a solução analítica duma equação diferencial num determinado intervalo, é o motivo principal para a introdução de métodos numéricos para aproximar essa solução. Estes métodos vão discretizar a equação diferencial nesse intervalo. As equações diferenciais envolvem a derivada das funções. São consideradas apenas equações diferenciais ordinárias (EDO), ie, com apenas uma variável independente. A ordem da maior derivada existente na equação diferencial define a sua ordem. Uma equação diferencial de ordem n precisa de n condições auxiliares para que a sua solução seja única. Estas condições auxiliares podem envolver o conhecimento da variável dependente no início do intervalo, tratando-se dum problema de condições iniciais. Quando as condições auxiliares dizem respeito a pontos na fronteira do intervalo, o problema diz-se de condições de fronteira. A sua distinção é muito importante já que os métodos para a sua resolução são completamente distintos. Vão resolver-se numericamente problemas com condições iniciais, utilizando o método de Runge-Kutta de 2a ordem para uma equação diferencial ordinária (A.33) e para um sistema de equações diferenciais ordinárias (A.33). Ainda nos problemas com condições iniciais, incluem-se problemas com equações diferenciais de ordem superior (n > 1). São também resolvidos problemas com condições de fronteira envolvendo diferenças finitas cen- trais (A.35), ascendentes (A.36) e descendentes (A.37). 149 151 q = hf(x2, y1 + p) = 1 × f(2, 83.041673) = 16.172054 y2 = 38.322923 + 1 2(44.718750 + 16.172054) = 68.768325 ≈ y(2) 3a etapa (i = 2): y3 = y2 + 1 2(p + q) p = hf(x2, y2) = 1 × [2 × 68.768325 × (92 − 68.768325 92 )] = 34.730508 q = hf(x3, y2 + p) = 1 × f(3, 103.498833) = −25.872083 y3 = 68.768325 + 1 2(34.730508 + (−25.872083)) = 73.197538 ≈ y(3) Ou seja, y1 = 38.322923 ≈ y(1), y2 = 68.768325 ≈ y(2), y3 = 73.197538 ≈ y(3). Ao fim de 3 horas, havia aproximadamente 73 alunos a terem conhecimento do boato. CONUM: EQUNIC Equações Diferenciais Ordinárias - Métodos de Passo Único Ordem = 2 Função f(x1,x2) = 2*x2*((92-x2)/92) Espaçamento entre Pontos = 1.000000 x0 = 0.000000 y0 = 10.000000 Limite Superior do Intervalo = 3.000000 Solução Numérica y0 = 10.000000 y1 = 38.322923 y2 = 68.768325 y3 = 73.197537 MATLAB: M-file: function dy=f7_1(x,y) dy=2*y*(92-y)/92; Comandos: >> [t,y]=ode45(’f7_1’,[0,1,2,3],[10]) t = y = 0 10.0000 1 43.6073 152 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 2 79.9871 3 90.1673 153 2. O progresso de uma epidemia de gripe numa população de N indivíduos é modelada pelo seguinte sistema de equações diferenciais:            dx dt = −βxy dy dt = βxy − αy dz dt = αy em que x é o número de pessoas susceptíveis de apanhar a gripe, y é o número de pessoas infectadas e z o número de pessoas imunes, incluindo todos os recuperados, no tempo t. Os parâmetros α e β são as taxas de recuperação e transmissão (por dia), respectivamente. Assume-se que a população é fixa, logo novos nascimentos são balanceados pelas mortes. Considere α = 0.05, β = 0.0002, x(0) = 980, y(0) = 20, z(0) = 0. a) Avalie a situação da população passados 10 dias do começo da epidemia, usando o método de Runge-Kutta de 2a ordem com h = 5. b) Interprete os resultados obtidos no fim de cada uma das etapas. Resolução: Mudança de variável: t → x, x → y1, y → y2 e z → y3. a) Trata-se de um sistema de 3 equações diferenciais de 1a ordem com condições iniciais. Discretizando a variável independente x com h = 5 obtém-se: x0 = 0, x1 = 5 e x2 = 10. As condições iniciais são as seguintes: y1,0 = y1(0) = 980, y2,0 = y2(0) = 20 e y3,0 = y3(0) = 0.          y′ 1 = −0.0002y1y2 = f1(x, y1, y2, y3) y′ 2 = 0.0002y1y2 − 0.05y2 = f2(x, y1, y2, y3) y′ 3 = 0.05y2 = f3(x, y1, y2, y3) 154 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA Este sistema será resolvido através do Método de Runge-Kutta de 2a ordem para sistemas (A.33). Em cada etapa, primeiro calculam-se os p′s e a seguir os q′s. 1a etapa (i = 0): p1 = hf1(x0, y1,0, y2,0, y3,0) = 5f1(0, 980, 20, 0) = 5 × (−0.0002 × 980 × 20) = −19.6 p2 = hf2(x0, y1,0, y2,0, y3,0) = 5f2(0, 980, 20, 0) = 5 × (0.0002 × 980 × 20 − 0.05 × 20) = 14.6 p3 = hf3(x0, y1,0, y2,0, y3,0) = 5f3(0, 980, 20, 0) = 5 × (0.05 × 20) = 5 q1 = hf1(x1, y1,0 + p1, y2,0 + p2, y3,0 + p3) = 5f1(5, 960.4, 34.6, 5) = −33.22984 q2 = hf2(x1, y1,0 + p1, y2,0 + p2, y3,0 + p3) = 5f2(5, 960.4, 34.6, 5) = 24.57984 q3 = hf3(x1, y1,0 + p1, y2,0 + p2, y3,0 + p3) = 5f3(5, 960.4, 34.6, 5) = 8.65 y1,1 = y1,0 + 1 2(p1 + q1) = 980 + 1 2(−19.6 + (−33.22984)) = 953.58508 ≈ y1(5) y2,1 = y2,0 + 1 2(p2 + q2) = 20 + 1 2(14.6 + 24.57984) = 39.58992 ≈ y2(5) y3,1 = y3,0 + 1 2(p3 + q3) = 0 + 1 2(5 + 8.65) = 6.825 ≈ y3(5) 2a etapa (i = 1): p1 = hf1(x1, y1,1, y2,1, y3,1) = 5f1(5, 953.58508, 39.58992, 6.825) = −37.752357 p2 = hf2(x1, y1,1, y2,1, y3,1) = 5f2(5, 953.58508, 39.58992, 6.825) = 27.854877 p3 = hf3(x1, y1,1, y2,1, y3,1) = 5f3(5, 953.58508, 39.58992, 6.825) = 9.89748 q1 = hf1(x2, y1,1 + p1, y2,1 + p2, y3,1 + p3) = 5f1(10, 915.832723, 67.444797, 16.72248) = −61.768152 q2 = hf2(x2, y1,1 + p1, y2,1 + p2, y3,1 + p3) = 5f2(10, 915.832723, 67.444797, 16.72248) = 44.906953 q3 = hf3(x2, y1,1 + p1, y2,1 + p2, y3,1 + p3) = 5f3(10, 915.832723, 67.444797, 16.72248) = 16.861199 y1,2 = y1,1 + 1 2(p1 + q1) = 953.58508 + 1 2(−37.752357 + (−61.768152)) = 903.824826 ≈ y1(10) y2,2 = y2,1 + 1 2(p2 + q2) = 39.58992 + 1 2(27.854877 + 44.906953) = 75.970835 ≈ y2(10) y3,2 = y3,1 + 1 2(p3 + q3) = 6.825 + 1 2(9.89748 + 16.861199) = 20.204340 ≈ y3(10) CONUM: SISDIF Equações Diferenciais Ordinárias - Sistemas de Equações Diferenciais Vector de Funções: -0.0002*x2*x3 ; 0.0002*x2*x3-0.05*x3 ; 0.05*x3 Valor Inicial: 0.000000 Vector Inicial dos y’s: 980.000000 ; 20.000000 ; 0.000000 Limite Superior do Intervalo: 10.000000 155 Espaçamento entre Pontos: 5.000000 Runge-Kutta Solução Numérica y1,0 = 980.000000 y2,0 = 20.000000 y3,0 = 0.000000 y1,1 = 953.585080 y2,1 = 39.589920 y3,1 = 6.825000 y1,2 = 903.824825 y2,2 = 75.970835 y3,2 = 20.204340 b) Ao longo do tempo, o número de pessoas susceptíveis de apanharem gripe (x - representado por y1) diminui (980 → 953.58508 → 903.824826), enquanto o número de infectados (y - representado por y2) aumenta (20 → 39.58992 → 75.970835). Ao mesmo tempo, o número de pessoas imunes (z - representado por y3) aumenta (0 → 6.825 → 20.204340). MATLAB: M-file: function dy=f7_2(x,y) dy=zeros(3,1); dy(1)=-0.0002*y(1)*y(2); dy(2)=0.0002*y(1)*y(2)-0.05*y(2); dy(3)=0.05*y(2); Comandos: >> [t,y]=ode45(’f7_2’,[0,5,10],[980,20,0]) t = 0 5 10 y = 980.0000 20.0000 0 951.6867 40.9841 7.3292 897.2762 80.6766 22.0472 156 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 3. Uma rolha de cortiça de comprimento L vai ser expulsa duma garrafa contendo um líquido em fermentação. As equações do movimento da rolha podem ser descritas pelas seguintes equações diferenciais:            dx dt = v dv dt =      g(1 + q)[(1 + x d)−γ + Rt 100 − 1 + qx L(1 + q)], x < L 0, x ≥ L em que g é a aceleração da gravidade (9.81 ms−2) q é o coeficiente de atrito da rolha x é o deslocamento da rolha no gargalo da garrafa t é o tempo d é o comprimento do gargalo da garrafa R é a razão percentual de aumento da pressão γ é a constante adiabática para o gas na garrafa (1.4) Considerando q = 20, L = 3.75cm, d = 5cm, R = 4, o sistema, após a substituição das constantes transforma-se em          dx dt = v dv dt =    206.01[(1 + 0.2x)−1.4 + 0.04t − 1 + 0.25397x], x < 3.75 0, x ≥ 3.75 As condições auxiliares do problema são x(0) = v(0) = 0. Enquanto que x < L a rolha mantém-se na garrafa, sendo expelida quando x = L. Considere o passo h = 0.75. Será que ao fim de duas etapas a rolha já saiu da garrafa? Qual a velocidade atingida na segunda etapa? Resolução: Mudança de variável: t → x, x → y1 e v → y2 157 Trata-se de um sistema de equações diferenciais de 1a ordem com condições iniciais (y1,0 = y1(0) = 0, y2,0 = y2(0) = 0). Discretizando em x, com h = 0.75: x0 = 0, x1 = 0.75 e x2 = 1.5. Substituindo as constantes, vem:          y′ 1 = f1(x, y1, y2) = y2 y′ 2 = f2(x, y1, y2) =    206.01[(1 + 0.2y1)−1.4 + 0.04x − 1 + 0.25397y1], y1 < 3.75 0, y1 ≥ 3.75 O sistema resolve-se pelo Método de Runge-Kutta de 2a ordem para sistemas (A.34): 1a etapa (i = 0): p1 = hf1(x0, y1,0, y2,0) = 0.75 × 0 = 0 p2 = hf2(x0, y1,0, y2,0) = 0.75 × 0 = 0 q1 = hf1(x1, y1,0 + p1, y2,0 + p2) = 0.75 × f1(0.75, 0, 0) = 0 q2 = hf2(x1, y1,0 + p1, y2,0 + p2) = 0.75 × f2(0.75, 0, 0) = 4.635225 y1,1 = y1,0 + 1 2(p1 + q1) = 0 + 1 2(0 + 0) = 0 ≈ y1(0.75) y2,1 = y2,0 + 1 2(p2 + q2) = 0 + 1 2(0 + 4.635225) = 2.317613 ≈ y2(0.75) 2a etapa (i = 1): p1 = hf1(x1, y1,1, y2,1) = 0.75 × f1(0.75, 0, 2.317613) = 1.738210 p2 = hf2(x1, y1,1, y2,1) = 0.75 × f2(0.75, 0, 2.317613) = 4.635225 q1 = hf1(x2, y1,1 + p1, y2,1 + p2) = 0.75 × f1(1.5, 1.738210, 6.952838) = 5.214629 q2 = hf2(x2, y1,1 + p1, y2,1 + p2) = 0.75 × f2(1.5, 1.738210, 6.952838) = 24.723309 y1,2 = y1,1 + 1 2(p1 + q1) = 0 + 1 2(1.738210 + 5.214629) = 3.476420 ≈ y1(1.5) y2,2 = y2,1 + 1 2(p2 + q2) = 2.317613 + 1 2(4.635225 + 24.723309) = 16.99688 ≈ y2(1.5) A variável y1 representa o espaço percorrido pela rolha. Como ao fim de 2 etapas y1,2 = 3.476420 é inferior ao comprimento da rolha, esta ainda não saiu da garrafa. y1,2 = x(1.5) = 3.476420 ≤ 3.75 (posição) A sua velocidade é y22 = v(1.5) = 16.99688. 158 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA MATLAB: M-file: function dy=f7_3(x,y) dy=zeros(2,1); dy(1)=y(2); if (y(1) < 3.75) dy(2)=206.01*[(1+0.2*y(1))^(-1.4)+0.04*x-1+0.25397*y(1)]; else dy(2)=0; end Comandos: >> [t,y]=ode45(’f7_3’,[0,0.75,1.5],[0,0]) t = 0 0.75 1.5 y= 0 0 0.535161 2.168723 7.931169 16.439906 159 4. Se y(t) representar a altitude de uma granada, então esta pode ser descrita pela seguinte equação diferencial: y′′(t) = −g+0.2y (g = 9.8m/s2) representa a aceleração da gravidade). Sabendo que a granada é lançada no instante t = 0, a partir do chão (y(0) = 0) e que, como a granada explode após 5 segundos, esta deverá estar a 40 m do chão. Será que, com uma velocidade inicial de 18 ms−1 (y′(0) = 18 ms−1), a gra- nada explode antes, depois ou na distância prevista? Utilize h = 2.5. Resolução: Mudança de variável: t → x Trata-se de uma equação diferencial de 2a ordem com condições iniciais. Vai surgir um sistema de 2 equações diferenciais de 1a ordem com condições iniciais, surgindo uma nova variável dependente. Fazendo y = y1 surge a nova variável dependente y2 = y′ 1    y′ 1 = y2 ≡ y′ 1 = f1(x, y1, y2) y′ 2 = −9.8 + 0.2y1 ≡ y′ 2 = f2(x, y1, y2) As condições iniciais para cada uma das variáveis dependentes são: y1,0 = y(0) = 0 e y2,0 = y2(0) = y′(0) = 18. O passo é de h = 2.5. Discretizando a variável independente x com h = 2.5: x0 = 0, x1 = 2.5 e x2 = 5. O sistema vai ser resolvido pelo Método de Runge-Kutta de 2a ordem para sistemas (A.34): 1a etapa (i = 0): p1 = hf1(x0, y1,0, y2,0) = 2.5f1(0, 0, 18) = 45 p2 = hf2(x0, y1,0, y2,0) = 2.5f2(0, 0, 18) = −24.5 q1 = hf1(x1, y1,0 + p1, y2,0 + p2) = 2.5f1(2.5, 45, −6.5) = −16.25 q2 = hf2(x1, y1,0 + p1, y2,0 + p2) = 2.5f2(2.5, 45, −6.5) = −2 160 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA y1,1 = y1,0 + 1 2(p1 + q1) = 0 + 1 2(45 + (−16.25)) = 14.375 ≈ y1(2.5) y2,1 = y2,0 + 1 2(p2 + q2) = 18 + 1 2(−24.5 + (−2)) = 4.75 ≈ y2(2.5) 2a etapa (i = 1): p1 = hf1(x1, y1,1, y2,1) = 2.5f1(2.5, 14.375, 4.75) = 11.875 p2 = hf2(x1, y1,1, y2,1) = 2.5f2(2.5, 14.375, 4.75) = −17.3125 q1 = hf1(x2, y1,1 + p1, y2,1 + p2) = 2.5f1(5, 26.25, −12.5625) = −31.40625 q2 = hf2(x2, y1,1 + p1, y2,1 + p2) = 2.5f2(5, 26.25, −12.5625) = −11.375 y1,2 = y1,1 + 1 2(p1 + q1) = 14.375 + 1 2(11.875 + (−31.40625)) = 4.609375 ≈ y1(5) y2,2 = y2,1 + 1 2(p2 + q2) = 4.75 + 1 2((−17.3125) + (−11.375)) = −9.59375 ≈ y2(5) Ao fim de 5 segundos, instante em que a granada explode, esta encontra-se aproxi- madamente a 4.609375 metros do ponto de partida, explodindo abaixo da altitude prevista (40 m). CONUM: SISDIF Equações Diferenciais Ordinárias - Sistemas de Equações Diferenciais Vector de Funções x3 -9.8+0.2*x2 Valor Inicial 0.000000 Vector Inicial dos y’s 0.000000 18.000000 Limite Superior do Intervalo 5.000000 Espaçamento entre Pontos 2.500000 Runge-Kutta Solução Numérica 161 y1,0 = 0.000000 y2,0 = 18.000000 y1,1 = 14.375000 y2,1 = 4.750000 y1,2 = 4.609375 y2,2 = -9.593750 Nota: Os resultados do CONUM são significativamente diferentes dos do MATLAB. No entanto, se o valor de h for diminuído (por exemplo h = 0.25), os resultados dos dois softwares aproximam-se. 162 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA MATLAB: M-file: function dy=f7_4(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=-9.8+0.2*y(1); Comandos: >> [t,y]=ode45(’f7_4’,[0,2.5,5],[0,18]) t = 0 2.5000 5.0000 y = 0 18.0000 21.0276 0.5390 3.2924 -16.1752 163 5. Um soldado pára-quedista cai do avião a uma altura de 600 metros. Após 5 segundos, o pára-quedas abre. A altura de queda do soldado pára-quedista como função do tempo, y(t), é dada por y′′ = −g + α(t) m , y(0) = 600m e y′(0) = 0m/s em que g = 9.81m/s2 é a aceleração da gravidade e m = 80 kg é peso do soldado pára-quedista. A resistência do ar α(t) é proporcional ao quadrado da velocidade, com diferentes constantes de proporcionalidade antes e depois da abertura do pára-quedas: α(t) =    K1y′(t)2, t < 5 s K2y′(t)2, t ≥ 5 s Considere K1 = 1/150, K2 = 4/150. A que altura o pára-quedas abre? (considere um espaçamento de 2.5 segundos). Resolução: Mudança de variável: t → x. Trata-se de uma equação diferencial de 2a ordem com condições iniciais, dando origem a um sistema de 2 equações de 1a ordem com condições iniciais. Vai surgir uma nova variável independente y = y1. As condições iniciais são: y1,0 = y1(0) = y(0) = 600 e y2,0 = y2(0) = y′ 1(0) = 0. Discretizando x, para h = 2.5, vem: x0 = 0, x1 = 2.5 e x2 = 5 (quando o pára- quedas abre). Substituindo as constantes, tem-se: Quando x < 5, utiliza-se a expressão:      y′ 1 = y2 = f1(x, y1, y2) y′ 2 = −9.81 + y2 2 80 × 150 = f2(x, y1, y2) 164 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA Quando x ≥ 5, usa-se:      y′ 1 = y2 = f1(x, y1, y2) y′ 2 = −9.81 + 4y2 2 80 × 150 = f2(x, y1, y2) O sistema de equações diferenciais vai ser resolvido usando o Método de Runge-Kutta de 2a ordem para sistemas (A.34): 1a etapa (i = 0): p1 = hf1(x0, y1,0, y2,0) = 2.5f1(0, 600, 0) = 0 p2 = hf2(x0, y1,0, y2,0) = 2.5f2(0, 600, 0) = −24.525 q1 = hf1(x1, y1,0 + p1, y2,0 + p2) = 2.5f1(2.5, 600, −24.525) = −61.3125 q2 = hf2(x1, y1,0 + p1, y2,0 + p2) = 2.5f2(2.5, 600, −24.525) = −24.399693 y1,1 = y1,0 + 1 2(p1 + q1) = 600 + 1 2(0 + (−61.3125)) = 569.34375 ≈ y1(2.5) y2,1 = y2,0 + 1 2(p2 + q2) = 0 + 1 2(−24.525 + (−24.399693)) = −24.462346 ≈ y2(2.5) 2a etapa (i = 1): NOTA: No cálculo dos q’s vai usar-se a expressão f2 para x ≥ 5 p1 = hf1(x1, y1,1, y2,1) = 2.5f1(2.5, 569.34375, −24.462346) = −61.155866 p2 = hf2(x1, y1,1, y2,1) = 2.5f2(2.5, 569.34375, −24.462346) = −24.400332 q1 = hf1(x2, y1,1 + p1, y2,1 + p2) = 2.5f1(5, 508.187884, −48.862678) = −122.156696 q2 = hf2(x2, y1,1 + p1, y2,1 + p2) = 2.5f2(5, 508.187884, −48.862678) = −22.535366 y1,2 = y1,1 + 1 2(p1 + q1) = 569.34375 + 1 2(−61.155866 + (−122.156696)) = 477.687469 ≈ y1(5) y2,2 = y2,1 + 1 2(p2 + q2) = −24.462347 + 1 2(−24.400332 + (−22.535366)) = −47.930195 ≈ y2(5) Nota: na 2a etapa, nos cálculos dos q’s utiliza-se a expressão para x ≥ 5. O pára-quedas abre aos aproximadamente 477.687469 metros de altitude. MATLAB: M-file: function dy=f7_5(x,y) dy=zeros(2,1); 165 dy(1)=y(2); if (x<5) aux=1/150*y(2)^2; else aux=4/150*y(2)^2; end dy(2)=-9.81+aux/80; Comandos: >> [t,y]=ode45(’f7_5’,[0,2.5,5],[600,0]) t = 0 2.5000 5.0000 y = 600.0000 0 569.3698 -24.4833 477.7904 -48.7108 166 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 6. Suponha um força f (f = 50) uniformemente distribuída ao longo do mastro de um barco. O vento origina a flexão do mastro. Para caracterizar este fenómeno e usando as leis da mecânica, tem-se a equação d2y dz2 = f 2 E I (L − z)2 em que E é o módulo de elasticidade (E = 1.5 × 108), L é a altura do mastro (L = 30) e I é o momento de inércia (I = 0.06). Em z = 0, y = 0 e dy dz = 0. Calcule numericamente a flexão no topo do mastro usando o método de Runge-Kutta de 2a ordem. Considere h = 15. Resolução: Mudança de variável: z → x. Trata-se de uma equação diferencial de 2a ordem com condições iniciais. Vai surgir um sistema de equações diferenciais de 1a ordem com condições iniciais: y1,0 = y1(0) = 0 e y2,0 = y2(0) = y′ 1(0) = 0. Discretizando a variável independente x com um passo de h = 15 obtém-se: x0 = 0, x1 = 15 e x2 = 30.    y′ 1 = y2 = f1(x, y1, y2) y′ 2 = 50 2 × 1.5 × 108 × 0.06(30 − x)2 = f2(x, y1, y2) O sistema será resolvido recorrendo ao Método de Runge-Kutta de 2a ordem para sistemas (A.34): 1a etapa (i = 0): 167 p1 = hf1(x0, y1,0, y2,0) = 15f1(0, 0, 0) = 0 p2 = hf2(x0, y1,0, y2,0) = 15f2(0, 0, 0) = 0.0375 q1 = hf1(x1, y1,0 + p1, y2,0 + p2) = 15f1(15, 0, 0.0375) = 0.5625 q2 = hf2(x1, y1,0 + p1, y2,0 + p2) = 15f2(15, 0, 0.0375) = 0.009375 y1,1 = y1,0 + 1 2(p1 + q1) = 0 + 1 2(0 + 0.5625) = 0.28125 ≈ y1(15) y2,1 = y2,0 + 1 2(p2 + q2) = 0 + 1 2(0.0375 + 0.009375) = 0.023438 ≈ y2(15) 2a etapa (i = 1): p1 = hf1(x1, y1,1, y2,1) = 15f1(15, 0.28125, 0.023438) = 0.35157 p2 = hf2(x1, y1,1, y2,1) = 15f2(15, 0.28125, 0.023438) = 0.009375 q1 = hf1(x2, y1,1 + p1, y2,1 + p2) = 15f1(30, 0.63282, 0.032813) = 0.492195 q2 = hf2(x2, y1,1 + p1, y2,1 + p2) = 15f2(30, 0.63282, 0.032813) = 0 y1,2 = y1,1 + 1 2(p1 + q1) = 0.28125 + 1 2(0.35157 + 0.492195) = 0.703133 ≈ y1(30) y2,2 = y2,1 + 1 2(p2 + q2) = 0.023438 + 1 2(0.009375 + 0) = 0.028126 ≈ y2(30) CONUM: SISDIF Equações Diferenciais Ordinárias - Sistemas de Equações Diferenciais Vector de Funções: x3 ; (50/(2*1.5e+08*0.06))*pot(30-x1,2) Valor Inicial: 0.000000 Vector Inicial dos y’s: 0.000000 ; 0.000000 Limite Superior do Intervalo: 30.000000 Espaçamento entre Pontos: 15.000000 Runge-Kutta Solução Numérica y1,0 = 0.000000 y2,0 = 0.000000 y1,1 = 0.281250 y2,1 = 0.023438 y1,2 = 0.703125 y2,2 = 0.028125 168 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA MATLAB: M-file: function dy=f7_6(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=50/(2*1.5*10^8*0.06)*(30-x)^2; Comandos: >> [t,y]=ode45(’f7_6’,[0,15,30],[0,0]) t = 0 15 30 y = 0 0 0.1992 0.0219 0.5625 0.0250 171 8. A equação de Duffing, descreve a dinâmica caótica de um circuito com um inductor não linear. A representação gráfica das variáveis de estado (y e dy dt ) ao longo do tempo (t), origina uma figura denominada mapa de Poincaré: d2y dt2 +ky+y3 = B cos t, Estime os valores das variáveis de estado, para k = 0.1, B = 12 e 0 ≤ t ≤ 0.1. Considere h = 0.05 e as condições y(0) = 0 e dy dt (0) = 4 usando o método de Runge- Kutta de 2a ordem. Resolução: Mudança de variável: t → x. Trata-se de uma equação diferencial de ordem 2 com condições iniciais. Vai dar origem a um sistema de 2 equações diferenciais de 1a ordem com condições iniciais. Discretizando x para h = 0.05 tem-se: x0 = 0, x1 = 0.05 e x2 = 0.1 (0 ≤ x ≤ 0.1).    y′ 1 = y2 = f1(x, y1, y2) y′ 2 = 12 cos(x) − 0.1y1 − y3 1 = f2(x, y1, y2) Usa-se o Método de Runge-Kutta de 2a ordem (A.34) com y1,0 = 0 e y2,0 = 4. 1a etapa (i = 0): p1 = hf1(x0, y1,0, y2,0) = 0.05f1(0, 0, 4) = 0.2 p2 = hf2(x0, y1,0, y2,0) = 0.05f2(0, 0, 4) = 0.6 q1 = hf1(x1, y1,0 + p1, y2,0 + p2) = 0.05f1(0.05, 0.2, 4.6) = 0.23 q2 = hf2(x1, y1,0 + p1, y2,0 + p2) = 0.05f2(0.05, 0.2, 4.6) = 0.597850 y1,1 = y1,0 + 1 2(p1 + q1) = 0 + 1 2(0.2 + 0.23) = 0.215 ≈ y1(0.05) y2,1 = y2,0 + 1 2(p2 + q2) = 4 + 1 2(0.6 + 0.597850) = 4.598925 ≈ y2(0.05) 2a etapa (i = 1): 172 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA p1 = hf1(x1, y1,1, y2,1) = 0.05f1(0.05, 0.215, 4.598925) = 0.229946 p2 = hf2(x1, y1,1, y2,1) = 0.05f2(0.05, 0.215, 4.598925) = 0.597678 q1 = hf1(x2, y1,1 + p1, y2,1 + p2) = 0.05f1(0.1, 0.444946, 5.196603) = 0.259830 q2 = hf2(x2, y1,1 + p1, y2,1 + p2) = 0.05f2(0.1, 0.444946, 5.196603) = 0.590373 y1,2 = y1,1 + 1 2(p1 + q1) = 0.215 + 1 2(0.229946 + 0.259830) = 0.459888 ≈ y1(0.1) y2,2 = y2,1 + 1 2(p2 + q2) = 4.598925 + 1 2(0.597678 + 0.590373) = 5.192951 ≈ y2(0.1) Os valores das variáveis de estado nos diferentes instantes são y1(0) = y1,0 = 0, y2(0) = y2,0 = 4, y1(0.05) ≈ y1,1 = 0.215, y2(0.05) ≈ y2,1 = 4.598925, y1(0.1) ≈ y1,2 = 0.459888 e y2(0.1) ≈ y2,2 = 5.192951. CONUM: SISDIF Equações Diferenciais Ordinárias - Sistemas de Equações Diferenciais Vector de Funções: x3 ; 12*cos(x1)-0.1*x2-pot(x2,3) Valor Inicial: 0.000000 Vector Inicial dos y’s: 0.000000 ; 4.000000 Limite Superior do Intervalo: 0.100000 Espaçamento entre Pontos: 0.050000 Runge-Kutta Solução Numérica y1,0 = 0.000000 y2,0 = 4.000000 y1,1 = 0.215000 y2,1 = 4.598925 y1,2 = 0.459888 y2,2 = 5.192951 MATLAB: M-file: function dy=f7_8(x,y) dy=zeros(2,1); dy(1)=y(2); dy(2)=12*cos(x)-0.1*y(1)-y(1)^3; Comandos: >> [t,y]=ode45(’f7_8’,[0,0.05,0.1],[0,4]) t = y = 0 0 4.0000 173 0.0500 0.2150 4.5991 0.1000 0.4598 5.1936 174 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 9. O potencial electrostático entre duas esferas de metal concêntricas de raios R1 e R2 (R1 < R2) é representado por u. O potencial da esfera interior é mantido constante com V1 volts e o potencial da esfera exterior é 0 volts. O potencial da região entre as duas esferas é governado pela equação de Laplace que, neste caso se reduz a d2u dr2 +2 r du dr = 0, para R1 ≤ r ≤ R2 com u(R1) = V1 e u(R2) = 0. Suponha que R1 = 2in, R2 = 4in e V1 = 110 volts. Calcule aproximações a u(2.5), u(3.0) e u(3.5). Resolução: Mudança de variável: r → x e u → y. Isolando o diferencial de maior ordem: y′′ + 2 xy′ = 0. Trata-se de uma equação diferencial de 2a ordem com condições de fronteira, pois as condições auxiliares dizem respeito aos valores de y na fronteira do intervalo. Discretizando a variável x (x ∈ [2, 4]) com h = 0.5: x0 = 2, x1 = 2.5, x2 = 3, x3 = 3.5, x4 = 4. Sabe-se que y(2) = 110 e y(4) = 0, que corresponde nesta discretização a y0 e y4, respectivamente. Existem 3 pontos interiores. Discretizando a equação, usando diferenças finitas centrais, nos pontos interiores (A.35), obtém-se: i = 1 1 h2 [y2 − 2y1 + y0] + 2 x1 ( 1 2h[y2 − y0]) = 0 1 0.52 [y2 − 2y1 + 110] + 2 2.5( 1 2 × 0.5[y2 − 110]) = 0 −8y1 + 4.8y2 = −352 175 i = 2 1 h2 [y3 − 2y2 + y1] + 2 x2 ( 1 2h[y3 − y1]) = 0 1 0.52 [y3 − 2y2 + y1] + 2 3.0( 1 2 × 0.5[y3 − y1]) = 0 3.333333y1 − 8y2 + 4.666667y3 = 0 i = 3 1 h2 [y4 − 2y3 + y2] + 2 x3 ( 1 2h[y4 − y2]) = 0 1 0.52 [0 − 2y3 + y2] + 2 3.5( 1 2 × 0.5[0 − y2]) = 0 3.428571y2 − 8y3 = 0 O sistema linear em y1, y2 e y3 é:      −8 4.8 0 | −352 3.333333 −8 4.666667 | 0 0 3.428571 −8 | 0      Resolvendo-se por EGPP, vem: y1 = 66 ≈ y(2.5), y2 = 36.666667 ≈ y(3) e y3 = 15.714286 ≈ y(3.5) CONUM: EQUFRO Equações Diferenciais Ordinárias - Problemas com Condições de Fronteira Coeficientes da Equação C: 0 ; 2/x1 ; 1 Coeficiente do Termo Independente: termo(x1) = 0 Coeficientes das Condições de Fronteira (limite inferior): 1 ; 0 Coeficientes das Condições de Fronteira (limite superior): 1 ; 0 Termos Independentes (inferior/superior): 110 ; 0 Valor Inicial: 2.000000 Limite Superior do Intervalo: 4.000000 Espaçamento entre Pontos: 0.500000 Solução Numérica y0 = 110.000000 y1 = 66.000000 y2 = 36.666667 y3 = 15.714286 y4 = 0.000000 176 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 10. A equação de Schrödinger da mecânica quântica pode ser escrita como: d2ψ(x) dx2 +[n(x)−βn]ψ(x) = 0 com as condições auxiliares dψ(−1) dx = −0.5, ψ(2) = 0 e x ∈ [−1, 2]. Considerando n(x) = n0 = 100 e βn = 95, construa e resolva o sistema linear de estimação de ψ(−1), ψ(0) e ψ(1). Resolução: Mudança de variável: ψ → y. Substituindo as constantes obtém-se: y′′(x) + 5y(x) = 0 com as condições auxiliares y′(−1) = −0.5 e y(2) = 0. Trata-se de uma equação diferencial de segunda ordem com condições de fronteira em [-1, 2], pois as condições auxiliares referem-se aos pontos de fronteira -1 e 2. O domínio será discretizado com h = 1, obtendo-se x0 = −1, x1 = 0, x2 = 1 e x3 = 2 (2 pontos interiores). Começa-se por substituir as derivadas da equação diferencial pelas respectivas dife- renças finitas centrais. A equação diferencial tem apenas uma segunda derivada, que será substituída por (A.35), ficando: 1 h2[yi+1 − 2yi + yi−1] + 5yi = 0 Substitui-se o i da equação anterior pelo índice dos pontos interiores (i = 1, 2): i = 1 → y0 + 3y1 + y2 = 0 i = 2 → y1 + 3y2 + y3 = 0 177 Obtém-se um sistema linear com 2 equações mas 4 incógnitas. Das condições de fronteira sabe-se que y3 = 0 pois y(x3) = y(2) = 0 e discretizando a equação de fronteira inferior com as diferenças finitas descendentes (A.37) obtém-se a equação: y′(−1) = −0.5 ≈ 1 1(y1 − y0) = −0.5 Colocando ordenadamente as 3 equações em sistema, vem:      −1 1 0 | −0.5 1 3 1 | 0 0 1 3 | 0      cuja solução, obtida por EGPP é: y0 = 0.363636 ≈ y(−1), y1 = −0.136364 ≈ y(0), y2 = 0.045455 ≈ y(1) Logo, ψ(−1) ≈ 0.363636, ψ(0) ≈ −0.136364 e ψ(1) ≈ 0.045455. CONUM: EQUFRO Equações Diferenciais Ordinárias - Problemas com Condições de Fronteira Coeficientes da Equação C: 5 ; 0 ; 1 Coeficiente do Termo Independente: termo(x1) = 0 Coeficientes das Condições de Fronteira (limite inferior): 0 ; 1 Coeficientes das Condições de Fronteira (limite superior): 1 ; 0 Termos Independentes (inferior/superior): -0.5 ; 0 Valor Inicial: -1.000000 Limite Superior do Intervalo: 2.000000 Espaçamento entre Pontos: 1.000000 Solução Numérica y0 = 0.363636 y1 = -0.136364 y2 = 0.045455 y3 = 0.000000 178 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 11. A equação que representa a perda de calor, T, por uma aleta é a seguinte: d2T dx2 −aT = −aTambiente com as seguintes condições: T0 = Tparede e dT(L) dx = 0. A segunda condição considera que a perda de calor na ponta da aleta é desprezável, pois o seu comprimento L é maior que a sua espessura. Para a = 20m−2o C−1, L = 0.3m, Tparede = 200o C e Tambiente = 20o C, estime a temperatura em 4 pontos igualmente distanciados da aleta. Resolução: Mudança de variável: T → y. Substituindo as constantes obtém-se: y′′ − 20y = −400, y(0) = 200, y′(0.3) = 0 Trata-se de uma equação diferencial de segunda ordem com condições de fronteira em [0, 0.3], pois as condições auxiliares referem-se aos pontos de fronteira 0 e 0.3. O domínio vai ser discretizado com h = 0.1, ficando x0 = 0, x1 = 0.1, x2 = 0.2 e x3 = 0.3. Começa-se por substituir as derivadas da equação diferencial pelas respectivas dife- renças finitas. Como a equação diferencial apenas tem uma segunda derivada, esta vai ser substituída por (A.35), ficando: 1 h2[yi+1 − 2yi + yi−1] − 20yi = −400 Substitui-se o i da equação anterior pelo índice dos pontos interiores (i = 1, 2): i = 1 → 100y0 − 220y1 + 100y2 = −400 179 i = 2 → 100y1 − 220y2 + 100y3 = −400 (sistema linear com 2 equações mas 4 incógnitas!) Das condições de fronteira tem-se que y0 = 200 que é substituído nas equações anteri- ores. Discretizando a equação fronteira superior com as diferenças finitas ascendentes (A.36) obtém-se a equação: y′(0.3) = 0 ≈ 1 0.1(y3 − y2) = 0 Colocando ordenadamente todas as equações em sistema, vem:      −220 100 0 | −20400 100 −220 100 | −400 0 −10 10 | 0      cuja solução, obtida por EGPP é:          y1 = 151.707317 ≈ y(0.1) y2 = 129.756098 ≈ y(0.2) y3 = 129.756098 ≈ y(0.3) CONUM: EQUFRO Equações Diferenciais Ordinárias - Problemas com Condições de Fronteira Coeficientes da Equação C: -20 ; 0 ; 1 Coeficiente do Termo Independente: termo(x1) = -400 Coeficientes das Condições de Fronteira (limite inferior): 1 ; 0 Coeficientes das Condições de Fronteira (limite superior): 0 ; 1 Termos Independentes (inferior/superior): 200 ; 0 Valor Inicial: 0.000000 Limite Superior do Intervalo: 0.300000 Espaçamento entre Pontos: 0.100000 Solução Numérica y0 = 200.000000 y1 = 151.707317 y2 = 129.756098 y3 = 129.756098 180 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 12. Em conforto térmico, é frequentemente analisada a perda de calor através de um determinado corpo, que é dada pela variação da temperatura T através da seguinte expressão: d2T dx −aT = −aText Num casaco de ski, de espessura 0.015 m e a = 0.13oC.m−2, com temperatura interior T0 = 28oC e temperatura exterior Text = −5oC, calcule a temperatura no interior do casaco em x = 0.005 e em x = 0.01. Resolução: Mudança de variável: T → y. Substituindo as constantes obtém-se: y′′ − 0.13y = 0.65 com as condições de fronteira y(0) = 28 e y(0.015) = −5. Trata-se de uma equa- ção diferencial de segunda ordem com condições de fronteira em [0, 0.015], pois as condições auxiliares referem-se aos pontos de fronteira 0 e 0.015. O domínio vai ser discretizado com h = 0.005, ficando x0 = 0, x1 = 0.005, x2 = 0.01 e x3 = 0.015. Começa-se por substituir as derivadas da equação diferencial pelas respectivas dife- renças finitas. Como a equação diferencial apenas tem uma segunda derivada, esta vai ser substituída por (A.35), ficando: 1 h2[yi+1 − 2yi + yi−1] − 0.13yi = 0.65 Substitui-se o i da equação anterior pelo índice dos pontos interiores (i = 1, 2): 40000(y2 − 2y1 + y0) − 0.13y1 = 0.65 ⇔ 40000y0 − 80000.13y1 + 40000y2 = 0.65 40000(y3 − 2y2 + y1) − 0.13y2 = 0.65 ⇔ 40000y1 − 80000.13y2 + 40000y3 = 0.65 Obtendo-se um sistema linear com 2 equações mas 4 incógnitas. 181 Das condições de fronteira tem-se que y0 = 28 e y3 = −5, que podem ser substi- tuídos nas equações anteriores, no entanto neste exercício optou-se por colocar estas condições também no sistema. Colocando todas as equações em sistema, vem:         1 0 0 0 | 28 40000 −80000.13 40000 0 | 0.65 0 40000 −80000.13 40000 | 0.65 0 0 0 1 | −5         cuja solução, obtida por EGPP é:                y0 = 28 = y(0) y1 = 16.999940 ≈ y(0.005) y2 = 5.999952 ≈ y(0.01) y3 = −5 = y(0.015) Assim, T1 ≈ T(0.005) ≈ 16.999940 e T2 ≈ T(0.015) ≈ 5.999952 oC. CONUM: EQUFRO Equações Diferenciais Ordinárias - Problemas com Condições de Fronteira Coeficientes da Equação C -0.13 ; 0 ; 1 Coeficiente do Termo Independente termo(x1) = 0.65 Coeficientes das Condições de Fronteira (limite inferior) 1 ; 0 Coeficientes das Condições de Fronteira (limite superior) 1 ; 0 Termos Independentes (inferior/superior) 28 ; -5 Valor Inicial 0.000000 Limite Superior do Intervalo 182 CAPÍTULO 8. DIFERENCIAÇÃO NUMÉRICA 0.015000 Espaçamento entre Pontos 0.005000 Matriz Tridiagonal do Sistema 1.000000 0.000000 0.000000 0.000000 40000.000000 -80000.130000 40000.000000 0.000000 0.000000 40000.000000 -80000.130000 40000.000000 0.000000 0.000000 -0.000000 1.000000 Termo Independente do Sistema 28.000000 0.650000 0.650000 -5.000000 Solução Numérica y0 = 28.000000 y1 = 16.999940 y2 = 5.999952 y3 = -5.000000 Bibliografia [1] CONUM, 2009. https://repositorium.sdum.uminho.pt/handle/1822/6017. [2] MATLAB, 2009. http://www.mathworks.com/. [3] E. Fernandes. COmputação NUMérica. Universidade do Minho, 1996. https:// repositorium.sdum.uminho.pt/handle/1822/5828. 183 Apêndice A Formulário • Fórmula fundamental dos erros δf ≤ Mx1δx1 + Mx2δx2 + ... + Mxnδxn (A.1) em que | ∂f ∂xi(ξ)| ≤ Mxi com ξ ∈ [x1 − δx1, x1 + δx1] × · · · × [xn − δxn, xn + δxn]. • Solução de uma equação não linear a) Método de Newton xk+1 = xk − f(xk) f ′(xk), k = 1, 2, ... (A.2) b) Método da Secante xk+1 = xk − (xk − xk−1)f(xk) f(xk) − f(xk−1) , k = 2, 3, ... (A.3) c) Critério de Paragem |xk+1 − xk| |xk+1| ≤ ǫ1 e |f(xk+1)| ≤ ǫ2 (A.4) • Métodos iterativos para sistemas de equações lineares 185 186 APÊNDICE A. FORMULÁRIO a) Método de Gauss-Seidel (D − L)x(k+1) = Ux(k) + b ou x(k+1) = CGSx(k) + (D − L)−1b (A.5) (A.6) CGS = (D − L)−1U b) Critério de Paragem ∥x(k+1) − x(k)∥ ∥x(k+1)∥ ≤ ǫ (A.7) • Sistemas de equações não lineares a) Método de Newton x(k+1) = x(k) + ∆x (A.8) J(x(k))∆x = −f(x(k)) b) Jacobiano J =      ∂f1(x1,x2,...,xn) ∂x1 . . . ∂f1(x1,x2,...,xn) ∂xn . . . . . . ∂fn(x1,x2,...,xn) ∂x1 . . . ∂fn(x1,x2,...,xn) ∂xn      (A.9) c) Critério de Paragem ∥∆x(k)∥ ∥x(k+1)∥ ≤ ǫ1 e ∥f(x(k+1))∥ ≤ ǫ2 (A.10) • Interpolação a) Tabela das diferenças divididas [xj, xj+1] = fj − fj+1 xj − xj+1 , j = 0, . . . , n − 1 (diferença dividida de ordem 1) (dd1) 190 APÊNDICE A. FORMULÁRIO y1,i+1 = y1,i + 1 2(p1 + q1), i = 0, 1... y2,i+1 = y2,i + 1 2(p2 + q2), i = 0, 1... . . . yn,i+1 = yn,i + 1 2(pn + qn), i = 0, 1.. p1 = h f1(xi, y1,i, y2,i, . . . , yn,i) q1 = h f1(xi+1, y1,i + p1, y2,i + p2, . . . , yn,i + pn) p2 = h f2(xi, y1,i, y2,i, . . . , yn,i) q2 = h f2(xi+1, y1,i + p1, y2,i + p2, . . . , yn,i + pn) ... ... pn = h fn(xi, y1,i, y2,i, . . . , yn,i) qn = h fn(xi+1, y1,i + p1, y2,i + p2, . . . , yn,i + pn) (A.34) c) Diferenças Finitas c1) Centrais y′(xi) ≈ 1 2h [y(xi+1) − y(xi−1)] y′′(xi) ≈ 1 h2 [y(xi+1) − 2y(xi) + y(xi−1)] y′′′(xi) ≈ 1 2h3 [y(xi+2) − 2y(xi+1) + 2y(xi−1) − y(xi−2)] yiv(xi) ≈ 1 h4 [y(xi+2) − 4y(xi+1) + 6y(xi) − 4y(xi−1) + y(xi−2)] (A.35) c2) Ascendentes y′(xi) ≈ 1 h [y(xi) − y(xi−1)] y′′(xi) ≈ 1 h2 [y(xi) − 2y(xi−1) + y(xi−2)] y′′′(xi) ≈ 1 h3 [y(xi) − 3y(xi−1) + 3y(xi−2) − y(xi−3)] (A.36)