Métodos Numéricos - norg.uminho.pt · Métodos Numéricos A. Ismael F. Vaz Departamento de Produção e Sistemas Escola de Engenharia Universidade do Minho [email protected] Mestrado

Métodos Numéricos

A. Ismael F. Vaz

Departamento de Produção e SistemasEscola de Engenharia

Universidade do [email protected]

Mestrado Integrado em Engenharia Mecânica

Ano lectivo 2007/2008

A. Ismael F. Vaz (UMinho) MN MIEMec. 2007/2008 1 / 203

Conteúdo

1 Introdução

2 Erros

3 Zeros de funções

4 Resolução de sistemas lineares

5 Resolução de sistemas não lineares

6 Interpolação polinomial

7 Interpolação segmentada - Splines

8 Mínimos quadrados lineares

9 Equações diferenciais - Condições iniciais

10 Equações diferenciais - Condições Fronteira

11 Integração numérica


Introdução

Conteúdo

1 Introdução

2 Erros











Introdução

Apresentação - Docente

Aulas teóricas

A. Ismael F. Vaz - [email protected]

www.norg.uminho.pt/aivaz

Aula teórico-prática e práticas

Senhorinha Teixeira – [email protected] Maria Rocha – [email protected]

Horário de atendimentoQuintas das 14h00 às 15h00.

Marcação por email.

As docentes das TPs e Ps terão o seu horário de atendimento.A. Ismael F. Vaz (UMinho) MN MIEMec. 2007/2008 4 / 203

Introdução


Aulas teóricas








Introdução


Aulas teóricas








Introdução

Apresentação - Disciplina

Página da disciplina;


6 fichas TPs a realizar ao longo do semestre (17 valores nas aulas Ts)e 6 fichas Ps a realizar ao longo do semestre (3 valores - 0.5 cada -nas aulas Ps).A classificação final é: Fichas TPs + Fichas Ps.

Não é obrigatória a presença nas aulas Ts, TPs ou Ps. Mas atençãoaos momentos de avaliação.A avaliação feitas nas aulas Ps é considerada avaliação laboratorial,pelo que é exigido o valor mínimo de 50% para ter frequência (poder ira exame).


Introdução







Introdução







Introdução







Introdução







Introdução

Material necessário e de apoio

Calculadora científica;

Folhas das fichas TPs;

Papel e caneta;Livro de Computação Numérica;www.norg.uminho.pt/emgpf

Software CoNum;www.norg.uminho.pt/emgpf

Formulário e acetatos disponíveis na página;www.norg.uminho.pt/aivaz


Introdução








Introdução








Introdução








Introdução








Introdução








Introdução

Programa detalhado

Dia Matéria28-Fev Apresentação da disciplina. Erros. Algarismos significativos. Fórmula

fundamental dos erros. Erros de truncatura.06-Mar Solução de equações não lineares. Método dos gráficos. Método da

secante e sua convergência. Método de Newton e sua convergência.Critérios de paragem.

13-Mar Sistemas de equações lineares. Eliminação de Gauss com pivotagemparcial. Avaliação sobre zeros de funções (2.5 valores).

27-Mar Métodos iterativos de Gauss-Seidel e Jacobi. Sistemas de equaçõesnão lineares. Método de Newton.

03-Abr Interpolação polinomial. Diferenças divididas. Fórmula interpoladorade Newton. Erro da fórmula interpoladora de Newton. Avaliaçãosobre sistemas lineares e não lineares (2.5 valores).


Introdução

Programa detalhado

Dia Matéria10-Abr Splines lineares e cúbicas.17-Abr Splines cúbicas. Revisões24-Abr Revisões. Avaliação sobre interpolação (3 valores).08-Mai Mínimos quadrados polinomiais e modelos lineares.30-Mai Equações diferenciais com condições iniciais de 1a. Sistemas de equa-

ções diferenciais de 1a ordem.05-Jun Equações diferenciais de ordem superior. Avaliação sobre mínimos

quadrados (3 valores).12-Jun Equações diferencias com condições fronteira. Integração numérica.

Fórmulas simples e compostas do Trapézio, Simpson e 3/8.19-Jun Revisões. Avaliação sobre equações diferenciais (3 valores).26-Jun Revisões. Avaliação sobre integração (3 valores).03-Jul Revisões.


Introdução

Motivação da disciplina

Presente em todos os cursos de engenharia (aplicações em todas asáreas da engenharia);

A disciplina de métodos numéricos dedica-se à resolução numérica deproblemas matemáticos. Com o desenvolvimento dos computadoresencontra-se direccionada para a implementação de algoritmos estáveis.


Introdução

Motivação da disciplina

Presente em todos os cursos de engenharia (aplicações em todas asáreas da engenharia);

A disciplina de métodos numéricos dedica-se à resolução numérica deproblemas matemáticos. Com o desenvolvimento dos computadoresencontra-se direccionada para a implementação de algoritmos estáveis.


Introdução

Controlo óptimo - Um exemplo

Problema de optimização do processo semi-contínuo de produção deEtanol.O problema de optimização é: (t0 = 0 e tf = 61.2 dias)

maxu(t)

J(tf ) ≡ x3(tf )x4(tf )

s.t.dx1

dt= g1x1 − u

x1

x4

dx2

dt= −10g1x1 + u

150− x2

x4

dx3

dt= g2x1 − u

x3

x4

dx4

dt= u

0 ≤ x4(tf ) ≤ 2000 ≤ u(t) ≤ 12∀t ∈ [t0, tf ]

com

g1 =(

0.4081 + x3/16

)(x2

0.22 + x2

)g2 =

(1

1 + x3/71.5

)(x2

0.44 + x2

)onde x1, x2 e x3 são as concentrações damassa celular, substrato e produto (g/L),e x4 é o volume (L). As condições iniciaissão:

x(t0) = (1, 150, 0, 10)T .


Introdução

Abordagem para a resolução

Grande exigência em termos numéricos;

Grande exigência em termos de programação;

Solução da equação diferencial com o CVODE (software em C);

Problemas codificados em AMPL (linguagem de modelação);

Algoritmo para optimização sem derivadas;


Introdução








Introdução








Introdução








Introdução








Erros

Conteúdo

1 Introdução

2 Erros











Erros

Formato de vírgula flutuante normalizado

fl(x) = ±0.d1d2...dk × 10e

onde, 0.d1d2 . . . dk corresponde à mantissa, e e é o expoente.flt(x) representa o valor de x em vírgula flutuante truncado efla(x) representa o valor de x em vírgula flutuante arredondado.

Exemplo

x =23

flt(x) = 0.66666× 100

fla(x) = 0.66667× 100


Erros

Formato de vírgula flutuante

(norma IEEE-754, 32 bits)

σ e + 64 d1 d2 d3 d4 d5 d6

1 bit 7 bits 4 bits 4 bits 4 bits 4 bits 4 bits 4 bits

Exemplox = 2490.125 = 9× 162 + 11× 161 + 10× 160 + 2× 16−1 =(9× 16−1 + 11× 16−2 + 10× 16−3 + 2× 16−4

)× 163

0 1000011 1001 1011 1010 0010 0000 0000σ e + 64 d1 d2 d3 d4 d5 d6


Erros

Exemplo de programação


Erros

Exemplo de programação


Erros

Erros

Seja x o valor exacto e x o seu valor aproximado, que será usado noscálculos

x− x é o erro absoluto (normalmente não se pode calcular, porque xé desconhecido);|x− x| ≤ δx é o limite superior do erro absoluto;

rx = |x−x||x| = δx

|x| ≈δx|x| é o erro relativo.

ExemploPediu-se a duas pessoas para contarem laranjas de dois cestos. A primeira contou980 laranjas num cesto de 1000 e a segunda contou 480 num cesto de 500.Apesar de cometerem o mesmo erro absoluto (δ1 = 20 laranjas e δ2 = 20 laranjas)a segunda cometeu um erro maior, visto que r1 = 20

1000 = 0.02 e r2 = 20500 = 0.04.


Erros

Erros



rx = |x−x||x| = δx



1000 = 0.02 e r2 = 20500 = 0.04.


Erros

Erros



rx = |x−x||x| = δx



1000 = 0.02 e r2 = 20500 = 0.04.


Erros

Erros



rx = |x−x||x| = δx



1000 = 0.02 e r2 = 20500 = 0.04.


Erros

Fórmula fundamental dos erros

Dados n valores aproximados, x1, . . . , xn, e os seus respectivos errosabsolutos é possível calcular um majorante para o erro absoluto cometidoquando se aplica uma função f , através da fórmula fundamental dos erros.

δf ≤Mx1δx1 + Mx2δx2 + ... + Mxnδxn

onde maxx∈I

∣∣∣ ∂f∂xi

∣∣∣ ≤Mxi , com I = Ix1 × · · · × Ixn eIxi = [xi − δxi , xi + δxi ]

rf ≤δf

|f(x1, . . . , xn)|


Erros

Exemplo

Cálculo dos limites do erro absoluto e relativo do cálculo da funçãof(x) = x1 − x2.

Temos que∣∣∣ ∂f∂x1

∣∣∣ ≤Mx1 = 1 e∣∣∣ ∂f∂x2

∣∣∣ ≤Mx2 = 1, logo

δf = δx1 + δx2

e

rf ≤δx1 + δx2

|x1 − x2|


Erros

Algarismos Significativos

Casa decimais são as casas (algarismos) à direita da vírgula.Os algarismos significativos são aqueles em que temos confiança do seuvalor.

Exemplos:0.1234567 tem 1 algarismo significativo se δ = 0.05, 2 se δ = 0.005 e 7 seδ = 0.00000005.

0.0000020 tem 7 casas decimais e 2 algarismos significativos(δ = 0.00000005).

Quando todas as casas decimais são significativas 0.2 é diferente de 0.20.


Zeros de funções

Conteúdo

1 Introdução

2 Erros











Zeros de funções

Forma geral do problema

Pretende-se determinar x∗ tal que

f(x) = 0

ExemploTemos x∗ = −0.567143290409784 como solução para

ex + x = 0

Nota: uma equação não linear pode não ter solução, ou ter mais do queuma.


Zeros de funções

Métodos iterativos

Uma sequência diz-se iterativa se é definida por uma função Findependente de k e dependente de um ou vários elementos anteriores a ele,

xk = F (xk−1, xk−2, . . . )

Aproximações iniciaisUm método que se baseie numa sequência iterativa com k − 1 elementosanteriores necessita também de k − 1 valores iniciais.

Exemploxk = xk−1 + xk−2

Partindo de x0 = 1 e x1 = 1 temos x2 = x1 + x0 = 2,x3 = x2 + x1 = 2 + 1 = 3, x4 = x3 + x2 = 3 + 2 = 5 gera uma sequênciacom os números de Fibonacci.


Zeros de funções

Convergência

Uma sequência iterativa diz-se convergente quando

limk→∞

xk = x∗

Convergência superlinear

limk→+∞

|x∗ − xk+1||x∗ − xk|1.618 = L ou lim

k→+∞

|x∗ − xk+1||x∗ − xk|

= 0

Convergência quadrática

limk→+∞

|x∗ − xk+1||x∗ − xk|2

= L


Zeros de funções

Critério de Paragem

A sequência de aproximações pode ser infinita. Como se pretende obteruma aproximação à solução implementa-se um critério de paragem.

Estimativa do erro relativo

dk =|xk+1 − xk||xk+1|

≤ ε1

Valor da função|f(xk+1)| ≤ ε2

Número máximo de iterações

k ≥ nmax


Zeros de funções

Método dos gráficos

Uma aproximação ao zero da função f(x) pode obter-se

pela intersecção do gráfico de f(x) com o eixo dos xx;se f(x) = g(x)− h(x) os zeros de f(x) são os pontos de intersecçãode g(x) com h(x).

O método dos gráficos é frequentemente usado para obtermos umaaproximação inicial para outros métodos mais precisos.


Zeros de funções

Exemplo

f(x) = ex + x g(x) = ex h(x) = −x

−1 −0.9 −0.8 −0.7 −0.6 −0.5 −0.4 −0.3 −0.2 −0.1 0−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

x

y

f(x)

h(x)

g(x)


Zeros de funções

Método da bissecção

Se f(xi)f(xs) < 0 então existe um número ímpar de raízes de f(x) nointervalo [xi, xs].

Aproxima-se da raiz calculando xk = xi+xs2 , k = 1, 2, . . .

Considera-se o intervalo

[xi, xk] se f(xi)f(xk) < 0 e faz-se xs ← xk

ou[xk, xs] se f(xk)f(xs) < 0 e faz-se xi ← xk


Zeros de funções

Interpretação gráfica (Bissecção)

f(x) = ex + x

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2−2

0

2

4

6

8

10

x

f(x)

xs

xi

xk

xs

xk+1


Zeros de funções

Método da secante

Método iterativo em que se fornece o x1 e x2 (a raiz não estánecessariamente no intervalo [x1, x2]). O próximo valor é calculado pelaseguinte fórmula (equação iterativa):

xk+1 = xk −(xk − xk−1)f(xk)f(xk)− f(xk−1)

, k = 2, 3, . . .

Zeros complexos: O método da secante também calcula zeros complexos,bastando para isso usar aritmética complexa.

Convergência: A convergência do método da Secante depende do valor deM2m ser pequeno. M é o max |f ′′(ξ)| e m é o min |f ′(η)|, onde ξ, η ∈ I.

εk+1 = − f ′′(ξ)2f ′(η)

εk−1εk


Zeros de funções

Interpretação gráfica (Secante)

f(x) = ex + x

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2−2

0

2

4

6

8

10

x

f(x)

xk xk−1 xk+1 xk+2


Zeros de funções

Método de Newton

Método iterativo em que se fornece o x0. O próximo valor é calculado pelaseguinte formula (equação iterativa):

xk+1 = xk −f(xk)f ′(xk)

, k = 1, 2, ...

Zeros complexos: O método de Newton também calcula zeroscomplexos, bastando para isso usar aritmética complexa.

Convergência: A convergência do método de Newton depende do valor deM2m ser pequeno. M é o max |f ′′(ξ)| e m é o min |f ′(η)|, onde ξ, η ∈ I.

εk+1 = − f ′′(ξ)2f ′(η)

ε2k


Zeros de funções

Interpretação gráfica (Newton)

f(x) = ex + x

−2 −1.5 −1 −0.5 0 0.5 1 1.5 2−2

0

2

4

6

8

10

x

f(x)

xk xk+1 xk+2


Zeros de funções

Principais propriedades

Ambos possuem convergência local. Superlinear no caso do métododa secante e quadrática no método de newton.

O método da secante não usa derivadas.

O método da secante e de Newton podem falhar quando odenominador da equação iterativa é próximo de zero, i.e., quandof(xk) ≈ f(xk−1) ou f ′(xk) ≈ 0.

O método da secante e de Newton não convergem necessariamentepara o zero mais próximo da aproximação inicial.


Zeros de funções

Exemplo ex + x = 0

Método de Newton com x0 = −0.5, ε1 = 0.5, ε2 = 0.1, nmax = 2.

Temos então que f(x) = ex + x e que f ′(x) = ex + 1.

1a iteração x0 = −0.5

f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)

O processo iterativo pára com x∗ ≈ x1 = −0.5665.

E se o ε1 fosse 0.1?A. Ismael F. Vaz (UMinho) MN MIEMec. 2007/2008 35 / 203

Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Zeros de funções

Exemplo ex + x = 0




f(−0.5) = e−0.5 − 0.5 = 0.1065 e f ′(−0.5) = 1.6065.

x1 = x0 − f(−0.5)f ′(−0.5) = −0.5− 0.1065

1.6065 = −0.5665

CP:

f(−0.5665) = 1.0082× 10−3 ≤ 0.1 (Verdadeiro)

|x1−x0||x1| = |−0.56665+0.5|

|−0.5665| = 0.1174 ≤ 0.5 (Verdadeiro)



Resolução de sistemas lineares

Conteúdo

1 Introdução

2 Erros













a11x1+ a12x2+ · · ·+ a1nxn = b1

a21x1+ a22x2+ · · ·+ a2nxn = b2

. . .an1x1+ an2x2+ · · ·+ annxn = bn

É um sistema com n equações lineares nas n incógnitas, x1, x2, . . . , xn. Osistema pode ser escrito na forma matricial Ax = b

a11 a12 . . . a1n

a21 a22 . . . a2n

. . .an1 an2 . . . ann

x1

x2

. . .xn

=

b1

b2

. . .bn



Exemplo

3 2 −29 7 −96 8 −8

x1

x2

x3

=

111

É um sistema linear de dimensão 3× 3. A matriz dos coeficientes

A =

3 2 −29 7 −96 8 −8

∈ R3×3 e o vector b = (1, 1, 1)T ∈ R3 é o vector dos

termos independentes.



Definições

A característica de uma matriz A, c(A), é o número máximo de linhasparalelas, ou colunas, linearmente independentes que existem namatriz.Para que um sistema seja possível e determinado temos de terc(A) = n. Caso contrário (c(A) < n) o sistema é indeterminado ouimpossível.À matrix (A|b) que se obtém ampliando A com a coluna do termoindependente b chama-se matriz ampliada do sistema.Triangular superior (inferior): É uma matriz em que os elementosabaixo (acima) da diagonal principal são zeros.Tridiagonal: Matriz em que aij = 0, se |i− j| ≥ 2, i, j = 1, . . . , n.Uma matriz com muitos elementos nulos diz-se esparsa.



Definições




Definições




Definições




Definições




Tipos de métodos

Métodos directos e estáveis. Métodos que calculam a solução exactado sistema ao fim de um número finito de operações elementares, casonão ocorram erros de arredondamento.

Matrizes dos coeficientes densas e de pequena dimensão.

Métodos iterativos. Métodos que definem uma sequência infinita deoperações, determinando uma sequência de aproximações, cujo limiteé a solução exacta do sistema.

Matrizes dos coeficientes esparsas e de grande dimensão.



Tipos de métodos

Métodos directos e estáveis. Métodos que calculam a solução exactado sistema ao fim de um número finito de operações elementares, casonão ocorram erros de arredondamento.

Matrizes dos coeficientes densas e de pequena dimensão.

Métodos iterativos. Métodos que definem uma sequência infinita deoperações, determinando uma sequência de aproximações, cujo limiteé a solução exacta do sistema.

Matrizes dos coeficientes esparsas e de grande dimensão.



Estabilidade numérica

Considere-se o seguinte sistema linear:{0.0001x1+ x2 = 1.0001

x1+ x2 = 2

cuja solução é x = (1, 1)T . Usando aritmética de três algarismossignificativos e considerando o multiplicador igual a− 1

0.100×10−3 = −0.100× 105, surge o sistema condensado{0.100× 10−3x1+ x2 = 0.100× 101

− 0.1× 105x2 = −0.1× 105

cuja solução é x = (0, 1)T !!!



Motivação - Continuação

Se nas mesmas condições usarmos a pivotagem parcial temos{x1+ x2 = 2

0.100× 10−3x1+ x2 = 0.100× 101

m = −0.100×10−3

1 = −0.100× 10−3

{x1+ x2 = 2

x2 = 0.100× 101

cuja solução é x = (1, 1)T .



Eliminação de Gauss com Pivotagem Parcial (EGPP)

Corresponde a eliminação de Gauss, mas em que a linha usada naeliminação dos elementos da coluna das linhas seguintes é o maior emmódulo.Exemplo: 3 2 −2

9 7 −96 8 −8

∣∣∣∣∣∣111

→ 9 7 −9

3 2 −26 8 −8

∣∣∣∣∣∣111

m21 = −39

m31 = −69 9 7 −9

0 −0.333333 10 3.333333 −2

∣∣∣∣∣∣1

0.6666670.333333

=



Eliminação de Gauss com Pivotagem Parcial (EGPP)

9 7 −90 3.333333 −20 −0.333333 1

∣∣∣∣∣∣1

0.3333330.666667

m32 = −−0.333333

3.333333 = 0.1=

9 7 −90 3.333333 −20 0 0.8

∣∣∣∣∣∣1

0.3333330.7



Substituição inversa

Quando a matriz é triangular superior pode-se determinar a soluçãodirectamente, através da substituição inversa.Exemplo 9 7 −9

0 3.333333 −20 0 0.8

∣∣∣∣∣∣1

0.3333330.7

vem que

x3 =0.70.8

= 0.875, x2 =0.333333− (−2)× 0.875

3.333333= 0.625

x1 =1− (−9)× 0.875− 7× 0.625

9= 0.5



Substituição directa

Quando a matriz é triangular inferior pode-se determinar a soluçãodirectamente, através da substituição directa.Exemplo 1 0 0

2 1 03 2 1

∣∣∣∣∣∣234

vem que

x1 =21

= 2, x2 =3− 2× 2

1= −1

x3 =4− 3× 2− 2× (−1)

1= 0



Decomposição LU

Da eliminação de Gauss com Pivotagem Parcial resulta

(A |I )→ (U |J )

Exemplo(1 2 1 02 1 0 1

)→(

2 1 0 11 2 1 0

)→(

2 1 0 10 1.5 1 −0.5

)



Determinantes de Matrizes

det(A) = (−1)sn∏

i=1

uii

onde uii corresponde aos elementos da diagonal da matriz U e s é onúmero de trocas de linhas para obter a matriz U .Exemplo

det

(1 22 1

)= (−1)1det

(2 10 1.5

)= (−1)1 × 2× 1.5 = −3



Cálculo da Inversa de Matrizes

A matriz inversa de A (A−1) verifica

AA−1 = I = A−1A.

O cálculo da matriz inversa reduz-se a resolução de n sistemas lineares daforma

Axj = ej , j = 1, . . . , n,

em que os vectores independentes ej são as colunas da matriz identidade.O vector solução xj corresponde à coluna j da matriz inversa.Na prática resolve-se os n sistemas em simultâneo, i.e., resolve-se aequação

(U |J )

por substituição inversa.



Cálculo da Inversa de Matrizes - Exemplo

(1 2 1 02 1 0 1

)→(

2 1 0 11 2 1 0

)→(

2 1 0 10 1.5 1 −0.5

)(

2 1 00 1.5 1

)→{

x11 = 0−1×0.66672 = −0.3334

x21 = 11.5 = 0.6667(

2 1 10 1.5 −0.5

)→{

x12 = 1−1×(−0.3333)2 = 0.6667

x22 = −0.51.5 = −0.3333

A inversa de(

1 22 1

)é(−0.3334 0.6667

0.6667 −0.3333

)



Métodos iterativos

Nos métodos iterativos a solução exacta só é obtida ao fim de umasequência infinita de operações.O processo parte de uma aproximação inicial para a solução do sistema eusa uma equação iterativa da forma

Mx(k+1) = Nx(k) + b, para k = 1, 2, . . .

Os métodos em que M e N não dependem de k dizem-se estacionários.Os métodos de Jacobi e Gauss-Seidel são métodos estacionário.



Método Iterativo Jacobi

D matriz dos elementos da diagonal principal, L matriz dos simétricos doselementos abaixo da diagonal principal e U matriz dos simétricos doselementos acima da diagonal principal.O método de Jacobi usa a partição de A em D − (L + U), i.e, M = D eN = L + UA equação iterativa fica

Dx(k+1) = (L + U)x(k) + b ou x(k+1) = D−1(L + U)x(k) + D−1b

A matriz iteração é

CJ = M−1N = D−1(L + U)



Método Iterativo Gauss-Seidel

M = D − L

N = U

A equação iterativa fica

Mx(k+1) = Nx(k) + b ou x(k+1) = M−1Nx(k) + M−1b

A matriz iteração é

CGS = M−1N = (D − L)−1U.



Critério de Paragem

Erro relativo na aproximação∥∥x(k+1) − x(k)∥∥∥∥x(k+1)

∥∥ < ε1

Resíduo ∥∥∥Ax(k+1) − b∥∥∥ < ε2



Convergência dos métodos iterativos

Condições suficientesA simétrica e definida positiva =⇒ GS exibe convergência global;A é estrita e diagonalmente dominante =⇒ J e GS exibemconvergência global;‖C‖p < 1, para qualquer normal p, =⇒ J e GS exibem convergênciaglobal;

C é a matriz iteração de Jacobi ou Gauss-Seidel.













Algumas definições

Uma matriz A diz-se simétrica se A = AT .Uma matriz é definida positiva se dT Ad > 0, ∀d 6= 0. É equivalente averificar que todos os determinante das sub-matrizes principais sãomaiores do que zero.Uma matriz A diz-se estrita e diagonalmente dominante se|aii| >

∑nj=1j 6=i

|aij |, i = 1, . . . , n





∑nj=1j 6=i

|aij |, i = 1, . . . , n





∑nj=1j 6=i

|aij |, i = 1, . . . , n



Exemplo - convergência de Gauss-Seidel

Considere-se a seguinte matriz dos coeficientes de um sistema linear

A =(

3 11 2

)Como a A = AT a matriz é simétrica.(

3 11 2

)det(|3|) = 3 > 0 det(A) = 3× 2− 1× 1 = 5 > 0

Logo A é simétrica e definida positiva e o método de Gauss-Seidelconverge.



Exemplo - convergência de Jacobi

Considere-se o seguinte sistema(1 2 13 1 1

)Como |1| ≯ |2| a matriz dos coeficientes não é estrita e diagonalmentedominante e nada se pode concluir acerca da convergência do método deJacobi. No entanto se trocarmos as linhas temos(

3 1 11 2 1

)e como |3| > |1| e |2| > |1| a matriz é estrita e diagonalmente dominante,logo o método de Jacobi converge.



Exemplo - convergência de Jacobi

Considere-se a seguinte matriz dos coeficientes de um sistema linear

A =(

3 23 1

)Como |3| > |2|, mas |1| ≯ |3| a matriz dos coeficientes não é estrita ediagonalmente dominante e nada se pode concluir acerca da convergênciado método de Jacobi.

D =(

3 00 1

)L =

(0 0−3 0

)U =

(0 −20 0

)



Continuação

CJ = D−1(L + U) =(

0.3333 00 1

)(0 −2−3 0

)=(

0 −0.6666−3 0

)Como

‖CJ‖∞ = max{|0|+ | − 0.6666|, | − 3|+ |0|} = 3 ≥ 1

e‖cJ‖1 = max{|0|+ | − 3|, | − 0.6666|+ |0|} = 3 ≥ 1

nada se pode concluir acerca da convergência do método de Jacobi.



Uma iteração do método de Gauss-Seidel

Considere-se o seguinte sistema linear

A =(

3 1 11 2 1

),

x(1) = (0, 0)T e ε1 = ε = 0.1

D =(

3 00 2

)L =

(0 0−1 0

)U =

(0 −10 0

)Equação iterativa é

(D − L)x(k+1) = Ux(k) + b



Continuação

1a iteração(3 01 2

)x(2) =

(0 −10 0

)(00

)+(

11

)=(

11

)(

3 0 11 2 1

)→

{x

(2)1 = 1

3 = 0.3333x

(2)2 = 1−1×0.3333

2 = 0.3334

C.P.

∥∥x(2) − x(1)∥∥∥∥x(2)

∥∥ =

∥∥∥∥( 0.33330.3334

)−(

00

)∥∥∥∥∞∥∥∥∥( 0.3333

0.3334

)∥∥∥∥∞

=0.33340.3334

= 1 ≮ 0.1

Como o critério não se verifica deve-se continuar com a próxima iteração.


Resolução de sistemas não lineares

Conteúdo

1 Introdução

2 Erros












Sistemas de equações não lineares

Forma geral do problemaf1(x1, x2, . . . , xn) = 0f2(x1, x2, . . . , xn) = 0. . .fn(x1, x2, . . . , xn) = 0

em que f = (f1, f2, . . . , fn)T é um vector de funções pelo menos uma vezcontinuamente diferenciáveis.Pretende-se determinar um x∗ = (x∗1, x

∗2, . . . , x

∗n)T tal que

f(x∗) = (0, 0, . . . , 0)T = 0.



Fórmula de Taylor a uma dimensão

Se f : R→ R for l + 1 vezes diferenciável temos que

f(x) =l∑

k=0

f (k)(a)k!

(x− a)k +f (l+1)(ξ)(l + 1)!

(x− a)l+1

com ξ ∈ [a, x] e a função definida em torno de a.

Exemplo: Valor da função em x(k+1) definido em torno de x(k).

f(x(k+1)) ≈ f(x(k)) + f ′(x(k))(x(k+1) − x(k))

ou seja, quando se pretende que f(x(k+1)) = 0 vem

x(k+1) = x(k) − f(x(k))f ′(x(k))

Eq. it. do método de Newton



Fórmula de Taylor para dimensão n

Se f : Rn → Rn temos que

f(x(k+1)) ≈ f(x(k))+∂f1(x(k))

∂x1

∂f1(x(k))∂x2

. . . ∂f1(x(k))∂xn

∂f2(x(k))∂x1

∂f2(x(k))∂x2

. . . ∂f2(x(k))∂xn

. . .∂fn(x(k))

∂x1

∂fn(x(k))∂x2

. . . ∂fn(x(k))∂xn

x(k+1)1 − x

(k)1

x(k+1)2 − x

(k)2

. . .

x(k+1)n − x

(k)n

e deduzindo a equação iterativa do método de Newton para sistemas deequações não lineares temos,

J(x(k))∆(k)x = −f(x(k)), com x(k+1) = x(k) + ∆(k)

x

em que J(x) é o Jacobiano da função.



Critério de paragem

‖x(k+1) − x(k)‖2‖x(k+1)‖2

=‖∆(k)

x ‖2‖x(k+1)‖2

≤ ε1

Se ‖x(k+1)‖2 é zero, ou próximo de zero, então o critério deve ser‖∆(k)

x ‖2 ≤ ε1

‖f(x(k+1))‖2 ≤ ε2

Número máximo de iterações.




‖x(k+1) − x(k)‖2‖x(k+1)‖2

=‖∆(k)

x ‖2‖x(k+1)‖2

≤ ε1


x ‖2 ≤ ε1

‖f(x(k+1))‖2 ≤ ε2





‖x(k+1) − x(k)‖2‖x(k+1)‖2

=‖∆(k)

x ‖2‖x(k+1)‖2

≤ ε1


x ‖2 ≤ ε1

‖f(x(k+1))‖2 ≤ ε2




Propriedades

Convergência local quadrática.

Determina a solução de um sistema linear numa única iteração.

Inconveniente do cálculo do Jacobiano. (Também existe um métododa secante para sistemas.)

O método falha quando o Jacobiano é singular (nova aproximaçãoinicial).

O método de Newton não converge necessariamente para a soluçãomais próxima da aproximação inicial.



Propriedades








Propriedades








Propriedades








Propriedades








Um exemplo

Considere-se o seguinte sistema não linear{3x2 + 2y2 = 354x2 − 3y2 = 24

cujo Jacobiano é J(x, y) =(

6x 4y8x −6y

)Temos

f(x, y) =(

3x2 + 2y2 − 354x2 − 3y2 − 24

)e a aproximação inicial é (x, y)(1) = (2.5, 2). Pretende-se determinar asolução com uma precisão de ε1 = ε2 = 10−1.



Continuação

1a iteração

J((x, y)(1)) = J(2.5, 2) =(

15 820 −12

)f((x, y)(1)) = f(2.5, 2) =

(−8.25−11

)

(15 8 8.2520 −12 11

)→(

20 −12 1115 8 8.25

)→(

20 −12 110 17 0

)→ ∆(1)

(x,y) =(

0.550

)e (x, y)(2) = (x, y)(1) + ∆(1)

(x,y) = (3.05, 2)T



Continuação

C.P. ∥∥∥f ((x, y)(2))∥∥∥

∞=∥∥∥∥( 0.9075

1.21

)∥∥∥∥∞

= 1.21 � ε2 = 0.1

Como o critério não se verifica faz-se uma nova iteração.


Interpolação polinomial

Conteúdo

1 Introdução

2 Erros












Motivação

Pretende-se determinar uma função aproximação que descreva o melhorpossível o comportamento de um conjunto de pontos (x0, f0), (x1, f1),. . . , (xm, fm).

Este conjunto de pontos pode ter sido obtido de:

observações de uma experiência (função desconhecida);uma função complexa cujo cálculo é difícil (função pode serconhecida).

A função aproximação server para:formular um modelo matemático que descreve o processo em causa;obter valores da função em pontos que são desconhecidos.

Problema: Como implementar a função sin(x) num microcontrolador?



Motivação








Motivação








Motivação








Motivação








Continuação

Pretende-se então, dado um conjunto de pontos (xi, fi), i = 1, . . . ,m,determinar uma função aproximação p(x) que melhor descreve ocomportamento dos dados, de acordo com uma certa medida.

No nosso caso vamos apenas considerar funções aproximação polinomiais,i.e., pn(x) é um polinómio interpolador de grau n.

Para construirmos o polinómio interpolador de Newton são necessárias asdiferenças divididas.



Diferenças divididas com espaçamento não constante

Considere-se uma função f(x) tabelada em m + 1 pontos x0, x1, . . . , xm

não igualmente espaçados.

Diferenças divididas de primeira ordem são

[xj , xj+1] =fj − fj+1

xj − xj+1j = 0, . . . ,m− 1

onde fj = f(xj).

A diferença dividida de primeira ordem corresponde ao declive da recta quepassa em (xj , fj) e (xj+1, fj+1).



Continuação

As diferenças divididas de segunda ordem são

[xj , xj+1, xj+2] =[xj , xj+1]− [xj+1, xj+2]

xj − xj+2, j = 0, . . . ,m− 2.

As diferenças divididas de ordem n são

[xj , xj+1, . . . , xj+n] =[xj , xj+1, . . . , xj+n−1]− [xj+1, xj+2, . . . , xj+n]

xj − xj+n

para j = 0, . . . ,m− n.



Tabela das diferenças divididas

x0 f0[x0, x1]

x1 f1 [x0, x1, x2][x1, x2] [x0, x1, x2, x3]

x2 f2 [x1, x2, x3][x2, x3] [x1, x2, x3, x4]

. . . . . . [x0, . . . , xm−1, xm]xm−2 fm−2 [xm−3, xm−2, xm−1]

[xm−2, xm−1] [xm−3, xm−2, xm−1, xm]xm−1 fm−1 [xm−2, xm−1, xm]

[xm−1, xm]xm fm



Exemplo

xi fi dd1 dd2 dd3 dd4 dd51 2

0.50003 3 0.1667

1.0000 −0.06674 4 −0.1667 0.0250

0.5000 0.0833 −0.00436 5 0.1667 −0.0139

1.0000 −0.01397 6 0.0833

1.333310 10



Propriedades das diferenças divididas

Simétrica nos argumentos, i.e., é independente da ordem dos argumentos;Exemplo

xi fi

6 51.0000

7 6 0.08331.3333

10 10

xi fi

7 61.3333

10 10 0.08331.25

6 5




Se fj = uj + vj para valores de xj , j = 0, . . . ,m então a tabela das DD def é igual à soma das tabelas das DD de u e v.Exemplo: f(x) = sin(x) + ex, u(x) = sin(x) e v(x) = ex.

xi ui

1 0.84150.0678

2 0.9093 −0.4180−0.7682

3 0.1411

xi vi

1 2.71834.6708

2 7.3891 4.012812.6964

3 20.0855

xi fi

1 3.55984.7386

2 8.2984 3.594811.9282

3 20.2266




A diferença dividida de cf(x), com c constante, é igual ao produto dadiferença dividida de f(x) por c.Exemplo: f(x) = sin(x), cf(x) = 2 sin(x)

xi fi

1 0.84150.0678

2 0.9093 −0.4180−0.7682

3 0.1411

xi 2fi

1 1.68300.1356

2 1.8186 −0.8360−1.5364

3 0.2822




As diferenças divididas de ordem n da função xn são iguais a 1 e as deordem r > n são nulas. Como consequência as diferenças divididas deordem n de um polinómio de ordem n são iguais e diferentes de zero.Exemplo: f(x) = x2 + 3x + 1

xi fi dd1 dd2 dd3 dd4−1 −1

20 1 1

4 01 5 1 0

6 02 11 1

83 19



Fórmula interpoladora de Newton

Das definições das diferenças divididas pode-se concluir que

f(x) = f0 + (x− x0)[x0, x] = f0 + (x− x0)f(x)− f0

x− x0

[x0, x] = [x0, x1] + (x− x1)[x0, x1, x][x0, x1, x] = [x0, x1, x2] + (x− x2)[x0, x1, x2, x]. . .

[x0, x1, . . . , xn−1, x] = [x0, x1, . . . , xn−1, xn] + (x− xn)[x0, x1, . . . , xn−1, xn, x]. . .

deduzindo-se a fórmula interpoladora de Newton

f(x) = f0 + (x− x0)[x0, x1] + (x− x0)(x− x1)[x0, x1, x2] + · · ·+(x− x0) . . . (x− xn−1)[x0, x1, . . . , xn] + · · ·



Polinómio interpolador de Newton

O polinómio interpolador de grau n obtém-se usado apenas n + 1 termosda fórmula interpoladora de Newton,

pn(x) = f0 + (x− x0)[x0, x1] + (x− x0)(x− x1)[x0, x1, x2] + · · ·+(x− x0) . . . (x− xn−1)[x0, x1, . . . , xn]

e temos queR(x) = f(x)− pn(x) = (x− x0)(x− x1) . . . (x− xn−1)(x− xn)f (n+1)(ξ)

(n+1)! ,ξ ∈ [x0, xn].



Diferenças divididas e derivadas

Da dedução da fórmula do erro de truncatura temos que

[x0, x1, . . . , xn] =f (n)(ξ)

n!, ξ ∈ [x0, xn],

ou seja, as diferenças divididas de primeira ordem são aproximações asprimeiras derivadas e a diferença dividida de ordem n é uma aproximação àderivada de ordem n da função sobre n!.



Determinação do polinómio interpolador

Em geral temos que m > n. Para construirmos o polinómio interpolador degrau n são necessários n + 1 pontos. A escolha dos pontos está relacionadacom o valor de x̄ para o qual se pretende obter uma estimativa da funçãof(x̄).

A escolha dos pontos deve obedecer as seguintes regras:os pontos xj e xj+1 em que xj < x̄ < xj+1 devem ser incluídos.os restantes, até formar os n + 1 necessários, são aqueles que estãomais próximos de x̄.



Determinação do polinómio interpolador

Em geral temos que m > n. Para construirmos o polinómio interpolador degrau n são necessários n + 1 pontos. A escolha dos pontos está relacionadacom o valor de x̄ para o qual se pretende obter uma estimativa da funçãof(x̄).

A escolha dos pontos deve obedecer as seguintes regras:os pontos xj e xj+1 em que xj < x̄ < xj+1 devem ser incluídos.os restantes, até formar os n + 1 necessários, são aqueles que estãomais próximos de x̄.



Exemplo

Considere-se a seguinte tabela de pontos

i 0 1 2 3 4 5xi 1 3 4 6 7 10fi 2 3 4 5 6 10

pretende-se obter uma estimativa de f(8) através usando uma interpolaçãoquadrática (polinómio de colocação de grau 2).Precisamos de 3 pontos para construir o polinómio de grau 2. x4 e x5

devem ser incluídos, uma vez que x4 < 8 < x5, e o restante será o x3.



Continuação

Construindo a tabela das diferenças divididas temos que

xi fi

6 51.0000

7 6 0.08331.3333

10 10

e p2(x) = 5 + (x− 6)× 1 + (x− 6)(x− 7)× 0.0833, ficandof(8) ≈ p2(8) = 5 + (8− 6)× 1 + (8− 6)(8− 7)× 0.0833 = 7.1666.Nota: x0 = 6, x1 = 7, e x2 = 10 para efeitos do cálculo do polinómio.Sendo o polinómio interpolador único qualquer combinação de pontosresulta no mesmo polinómio.



Continuação

xi fi

7 61.3333

10 10 0.08331.25

6 5

e p2(x) = 6 + (x− 7)× 1.3333 + (x− 7)(x− 10)× 0.0833, ficandof(8) ≈ p2(8) = 6 + (8− 7)× 1.333 + (8− 7)(8− 10)× 0.0833 = 7.1667.



Cálculo do majorante do erro absoluto

ET (x) =

∣∣∣∣∣(x− x0)(x− x1)(x− x2)f (3)(ξ)

3!

∣∣∣∣∣ , ξ ∈ [6, 10]

Como necessitámos de uma estimativa para f (3)(ξ) e a função f(x) édesconhecida vamos usar as diferenças divididas para a obter.

Com os pontos usados na construção do polinómio não é possível obteruma diferença dividida de ordem 3 e por isso vamos acrescentar mais umponto à tabela anterior.

É indiferente inserir o ponto no início ou no final da tabela, uma vez que ovalor de dd3 será o mesmo.



Continuação

xi fi

4 40.5000

6 5 0.16671.0000 −0.0139

7 6 0.08331.3333

10 10

ET (8) = |(8− 6)(8− 7)(8− 10)× (−0.0139)| = 0.0556


Interpolação segmentada - Splines

Conteúdo

1 Introdução

2 Erros












Motivação

Pretende-se aproximar uma função f(x) conhecida em alguns pontosdo seu domínio

por ser desconhecida a sua expressão analítica (experiêncialaboratorial);porque a sua expressão analítica é demasiado complexa.

Dará bons resultados aproximar uma função conhecida em 11 pontospor um polinómio de grau 10?



Motivação






Motivação






Motivação






Exemplo

Dada a função f(x) conhecida para os seguintes pontos (i = 0, . . . , 10)

xi 1 2 3 4 5 6 7 8 9 10 11fi 0.1 0.2 0.3 0.4 0.45 0.44 0.4 0.45 0.4 0.4 0.35

0 2 4 6 8 10 120

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

x

p 10(x

)



Interpolação segmentada

A ideia subjacente à interpolação segmentada consiste na construçãode segmentos entre pontos (em oposição à construção de um únicopolinómio).A facilidade de usar polinómios fazem com que sejam candidatosnaturais para a construção dos segmentos:

quando o segmento é um polinómio de grau 1 estamos na presença deuma Spline linear;quando o segmento é um polinómio de grau 3 estamos na presença deuma Spline cúbica.

O grau do polinómio usado na construção do segmento não estárelacionado com o número de pontos.Para um conjunto de n + 1 pontos temos n segmentos.

































Exemplo

0 2 4 6 8 10 120

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

x

p 10(x

), s

1(x),

s3(x

)

s110(x)

s31(x)

s32(x)

s33(x)

s19(x)

s18(x)

p10

(x)



Splines lineares

Considere-se dois pontos consecutivos (xi−1, xi) na construção de umsegmento da Spline linear. Podemos usar, por exemplo, o polinómiointerpolador de Newton para deduzir a equação da Spline linear.

Equação da Spline linear

xi−1 fi−1fi−fi−1

xi−xi−1si1(x) = fi−1 + fi−fi−1

xi−xi−1(x− xi−1)

xi fi

Interpolação nos segmentos

O segmento si1(x) apenas está definido para x ∈ [xi−1, xi], i = 1, . . . , n.



Splines lineares

Considere-se dois pontos consecutivos (xi−1, xi) na construção de umsegmento da Spline linear. Podemos usar, por exemplo, o polinómiointerpolador de Newton para deduzir a equação da Spline linear.

Equação da Spline linear

xi−1 fi−1fi−fi−1

xi−xi−1si1(x) = fi−1 + fi−fi−1

xi−xi−1(x− xi−1)

xi fi

Interpolação nos segmentos

O segmento si1(x) apenas está definido para x ∈ [xi−1, xi], i = 1, . . . , n.



Propriedades

Função segmentadaA Spline linear é uma função definida por segmentos.

s1(x) =

s11(x) se x ∈ [x0, x1]

s21(x) se x ∈ [x1, x2]· · ·sn1 (x) se x ∈ [xn−1, xn]

Mas na grande maioria das vezes só se calcula o segmento necessário.

Simplicidade vs diferenciabilidadeSão segmentos simples de determinar, mas a função final, no caso geral,não é diferenciável (nos nós intermédios).



Propriedades

Função segmentadaA Spline linear é uma função definida por segmentos.

s1(x) =

s11(x) se x ∈ [x0, x1]

s21(x) se x ∈ [x1, x2]· · ·sn1 (x) se x ∈ [xn−1, xn]

Mas na grande maioria das vezes só se calcula o segmento necessário.

Simplicidade vs diferenciabilidadeSão segmentos simples de determinar, mas a função final, no caso geral,não é diferenciável (nos nós intermédios).



Erro de truncatura

Erro de truncatura - spline linear

|f(x)− s1(x)| ≤ 18h2M2

Onde h = maxi(xi − xi−1) e M2 é um majorante do valor absoluto def ′′(x) em [x0, xn] (maxξ∈[x0,xn] |f ′′(ξ)| ≤M2).

Majorante de f ′′(x)

Quando não é possível obter o majorante de outra forma deve-se usar asdiferenças divididas para obter uma estimativa (atenção ao 2!).



Erro de truncatura

Erro de truncatura - spline linear

|f(x)− s1(x)| ≤ 18h2M2

Onde h = maxi(xi − xi−1) e M2 é um majorante do valor absoluto def ′′(x) em [x0, xn] (maxξ∈[x0,xn] |f ′′(ξ)| ≤M2).

Majorante de f ′′(x)

Quando não é possível obter o majorante de outra forma deve-se usar asdiferenças divididas para obter uma estimativa (atenção ao 2!).



Splines cúbicas

Considere-se um conjunto de pontos x0 ≤ x1 ≤ · · · ≤ xn.

DefiniçãoUm segmento cúbico (polinómio de grau 3) entre dois pontos está bemdefinido?

Resposta

Claro que não (existe um infinidade de polinómios de grau 3 queinterpolam os dois pontos).

Condições necessárias para definir a spline cúbica

Serão necessárias condições adicionais (de continuidade da spline) para quese consiga deduzir a(s) sua(s) equação(ões).



Splines cúbicas



Resposta






Splines cúbicas



Resposta






Dedução da equação

si3(x) é um polinómio de grau 3, logo a sua segunda derivada é uma recta.

Sejam M(x0) = M0, M(x1) = M1, . . . , M(xn) = Mn as segundasderivadas da spline cúbica nos nós x0, x1, . . . , xn, respectivamente.

Então temos que

(si3(x))′′

(si3(x))′′ =

x− xi

xi−1 − xiMi−1 +

x− xi−1

xi − xi−1Mi

A spline cúbica pode também ser usada para aproximar as segundasderivadas.



Equação da spline

Integrando duas vezes (si3(x))′′ obtém-se

si3(x), i = 1, . . . , n

si3(x) =

Mi−1

6(xi − xi−1)(xi − x)3 +

Mi

6(xi − xi−1)(x− xi−1)3+[

fi−1

(xi − xi−1)− Mi−1(xi − xi−1)

6

](xi − x)+[

fi

(xi − xi−1)− Mi(xi − xi−1)

6

](x− xi−1)

Impondo as condições de continuidade da spline e o facto de que cadasegmento é um polinómio interpolador (por exemplo si

3(xi) = fi).



Condições adicionais

Os valores de Mi, i = 0, . . . , n são, usualmente, incógnitas (segundasderivadas da função f(x) nos nós).

Continuidade da primeira derivada da spline

Resulta que (si3(xi))′ = (si+1

3 (xi))′, i = 1, . . . , n− 1, i.e.,

(xi − xi−1)Mi−1 + 2(xi+1 − xi−1)Mi + (xi+1 − xi)Mi+1 =6

xi+1 − xi(fi+1 − fi)−

6xi − xi−1

(fi − fi−1)

para o nó i (nós interiores).

Temos então n− 1 equações lineares para determinar n + 1 incógnitas, oque torna o sistema linear indeterminado (infinidade de soluções).



Tipos de splines cúbicas

Natural

M0 = Mn = 0, i.e., (s13(x0))′′ = (sn

3 (xn))′′ = 0, i.e., curvatura nula nasextremidades

Completa

Impondo que (s13(x0))′ = f ′(x0) e (xn

3 (xn))′ = f ′(xn) resulta em maisduas equações:

2(x1 − x0)M0 + (x1 − x0)M1 =6

x1 − x0(f1 − f0)− 6f ′(x0)

2(xn − xn−1)Mn + (xn − xn−1)Mn−1 = 6f ′(xn)− 6xn − xn−1

(fn − fn−1)

LivreNa prática desde que sejam dadas mais duas condições é possíveldeterminar a spline de forma única.




Natural

M0 = Mn = 0, i.e., (s13(x0))′′ = (sn


Completa



2(x1 − x0)M0 + (x1 − x0)M1 =6

x1 − x0(f1 − f0)− 6f ′(x0)


(fn − fn−1)





Natural

M0 = Mn = 0, i.e., (s13(x0))′′ = (sn


Completa



2(x1 − x0)M0 + (x1 − x0)M1 =6

x1 − x0(f1 − f0)− 6f ′(x0)


(fn − fn−1)




Spline cúbica completa

Caso em que derivadas são desconhecidasNo caso em que se pretende uma spline cúbica completa, mas as derivadasnos extremos são desconhecidas, estas podem ser estimadas por diferençasdivididas usando dois pontos auxiliares a e b.

f ′(x0) =f(a)− f0

a− x0f ′(xn) =

f(b)− fn

b− xn

Neste caso os pontos usados na obtenção da estimativa das derivadasdevem ser removidas da construção da spline (usualmente considerámosa = x1 e b = xn−1).

Se conhecer f(x)

Se f(x) for conhecida deve-se usar a derivada para obter os valores(exactos) de f ′0 e f ′n.



Spline cúbica completa

Caso em que derivadas são desconhecidasNo caso em que se pretende uma spline cúbica completa, mas as derivadasnos extremos são desconhecidas, estas podem ser estimadas por diferençasdivididas usando dois pontos auxiliares a e b.

f ′(x0) =f(a)− f0

a− x0f ′(xn) =

f(b)− fn

b− xn

Neste caso os pontos usados na obtenção da estimativa das derivadasdevem ser removidas da construção da spline (usualmente considerámosa = x1 e b = xn−1).

Se conhecer f(x)

Se f(x) for conhecida deve-se usar a derivada para obter os valores(exactos) de f ′0 e f ′n.



Spline

Spline cúbicaCada segmento i só pode ser usado para interpolar a função no intervalo[xi−1, xi], uma vez que a spline é uma função definida por segmentos:

s3(x) =

s13(x) se x ∈ [x0, x1]

s23(x) se x ∈ [x1, x2]· · ·

sn3 (x) se x ∈ [xn−1, xn]

Que segmento?Na prática determinamos apenas o segmento que possa interessar paraobter uma estimativa da função.



Spline

Spline cúbicaCada segmento i só pode ser usado para interpolar a função no intervalo[xi−1, xi], uma vez que a spline é uma função definida por segmentos:

s3(x) =

s13(x) se x ∈ [x0, x1]

s23(x) se x ∈ [x1, x2]· · ·

sn3 (x) se x ∈ [xn−1, xn]

Que segmento?Na prática determinamos apenas o segmento que possa interessar paraobter uma estimativa da função.



Erro de truncatura

Erro de truncatura da spline

|f(x)− s3(x)| ≤ 5384

h4M4

em que h é o maior espaçamento entre pontos usados na construção daspline e M4 e um majorante do valor absoluto da quarta derivada de f(x)no intervalo [x0, xn], i.e,

maxξ∈[x0,xn]

|f (iv)(ξ)| ≤M4

Erro de truncatura da derivada da spline

|f ′(x)− s′3(x)| ≤ 124

h3M4



Erro de truncatura

Erro de truncatura da spline

|f(x)− s3(x)| ≤ 5384

h4M4

em que h é o maior espaçamento entre pontos usados na construção daspline e M4 e um majorante do valor absoluto da quarta derivada de f(x)no intervalo [x0, xn], i.e,

maxξ∈[x0,xn]

|f (iv)(ξ)| ≤M4

Erro de truncatura da derivada da spline

|f ′(x)− s′3(x)| ≤ 124

h3M4



Árvore de decisão

Spline cúbica

Natural

Usando todos ospontos da tabela,considerar M0 =Mn = 0.

Completa

Deriv. desc. Deriv. conhec.f ′(x0) =?, f ′(xn) =? f ′(x0), f ′(xn)

Estimar as derivadaspor diferenças dividi-das, removendo doispontos na tabela

Usar as derivadas nasduas equações adicio-nais

Resolver um sistema de equações lineares por EGPP para obter os M ’s.

Determinar segmento e calcular aproximação.



Exemplo

EnunciadoAo efectuar observações astronómicas medindo as variações na magnitudeaparente, M , de uma estrela variável chamada variável Cepheid, ao longode um período de tempo, t, foram obtidos os seguintes valores:

tempo (t) 0.0 0.3 0.5 0.6 0.8Magnitude aparente (M) 0.302 0.106 0.240 0.579 0.468

Determine um valor aproximado da magnitude aparente da variável Cepheidno instante t = 0.4, utilizando uma spline cúbica completa.

Como se trata de uma spline completa em que se desconhece as derivadasem t = 0 e t = 0.8, considerámos a = 0.3 e b = 0.6 como pontos auxiliaresno cálculo da aproximação a f ′0 e f ′n (mudança de notação para não darconfusão com as incógnitas da spline).



Resolução

Temos então que:

f ′0 =0.106− 0.302

0.3− 0= −0.6533 f ′2 =

0.579− 0.4680.6− 0.8

= −0.555

e ficámos comi 0 1 2ti 0.0 0.5 0.8fi 0.302 0.240 0.468

Só existe uma equação para o único nó interior (i = 1).

(t1 − t0)M0 + 2(t2 − t0)M1 + (t2 − t1)M2 =6

t2 − t1(f2 − f1)−

6t1 − t0

(f1 − f0)

0.5M0 + 1.6M1 + 0.3M2 =6(0.468− 0.240)

0.3− 6(0.240− 0.302)

0.5= 5.304



Resolução – Cont.

Como a spline é completa temos as duas condições adicionais:

2(t1 − t0)M0 + (t1 − t0)M1 =6

t1 − t0(f1 − f0)− 6f ′0

2(t2 − t1)M2 + (t2 − t1)M1 = 6f ′2 −6

t2 − t1(f2 − f1)

Substituindo:

M0 + 0.5M1 =6(0.240− 0.302)

0.5− 6× (−0.6533) = 3.1758

0.6M2 + 0.3M1 = 6× (−0.555)− 6(0.468− 0.240)0.3

= −7.89

Resultando no seguinte sistema linear: 1.0 0.5 0.0 3.17580.5 1.6 0.3 5.30400.0 0.3 0.6 −7.89

→

M0 = −0.0163M1 = 6.3843M2 = −16.3421



Segmento

Como t = 0.4 se encontra no primeiro segmento da spline temos:

s13(t) =

M0

6(t1 − t0)(t1 − t)3 +

M1

6(t1 − t0)(t− t0)3+[

f0

(t1 − t0)− M0(t1 − t0)

6

](t1 − t)+[

f1

(t1 − t0)− M1(t1 − t0)

6

](t− t0)

s13(t) =

−0.01633

(0.5− t)3 +6.3843

3t3+[

0.3020.5

− −0.0163× 0.56

](0.5− t)+[

0.2400.5

− 6.3843× 0.56

]t



Aproximação

s13(t) = −0.0005(0.5− t)3 + 2.1281t3 + 0.6054(0.5− t)− 0.0520t

Não desenvolver o polinómio

s13(0.4) = −0.0005(0.5− 0.4)3 + 2.1281× 0.43

+ 0.6054(0.5− 0.4)− 0.0520× 0.4 = 0.1759



Erro de truncatura

Para o cálculo do erro de truncatura da spline temos calcular

h = maxi

(ti − ti−1) = max(0.5− 0.0, 0.8− 0.5) = max(0.5, 0.3) = 0.5

Uma vez que a função é desconhecida (majorante da quarta derivada)temos que usar as diferenças divididas.

ti fi(Mi) dd1 dd2 dd3 dd40.0 0.302

−0.65330.3 0.106 2.6467

0.6700 10.70000.5 0.240 9.0667 -68.9167

3.3900 −44.43330.6 0.579 −13.1500

−0.55500.8 0.468



Erro de truncatura

Logo temos que M4 = | − 68.9167× 4!| = 68.9167× 4!, resultando em

(Se existissem mais diferenças divididas de quarta hora escolher-se-ia a demaior valor absoluto e se não existirem pontos suficientes não é possívelcalcular o erro de truncatura.)

|s3(x)− f(x)| ≤ 5384× 0.54 × 68.9167× 4! = 1.3460


Mínimos quadrados lineares

Conteúdo

1 Introdução

2 Erros













MotivaçãoOs mínimos quadrados são usados quando os dados obtidos para umadeterminada função f(x) estão afectados de erros (ou ruídos). Neste casonão faz muito sentido usar uma interpolação, mas sim construir umafunção que reflicta, na generalidade, o comportamento dos dados.

ExemploSuponhamos que se pretende estimar o consumo de um automóvel emfunção da velocidade. Através da realização de várias experiênciaschegou-se aos seguintes valores:

xi(velocidade - km/h) 20 40 50 60 60 70 80fi(consumo l/100km) 5.5 5.6 5.7 5.9 5.85 6.1 7.5



Cont.

Seria correcto usar um polinómio interpolador para modelar a função f(x)?

20 30 40 50 60 70 805.4

5.6

5.8

6

6.2

6.4

6.6

6.8

7

7.2

7.4

Velocidade (km/h)

Con

sum

o (l/

100k

m)

O mais correcto seria determinar um polinómio de grau 2 (ou c1 + c2ex?)

que melhor se aproxima-se dos dados da tabela.



Forma geral do problema - Caso discreto

Dado um conjunto de valores (xi, fi), i = 1, . . . ,m, pretende-se determinarum modelo M(x) que aproxima o melhor possível a função dada, ou seja,

minimizar〈f −M(x), f −M(x)〉

em que

〈g(x), h(x)〉 =m∑

i=1

ω(xi)g(xi)h(xi).

Pretende-se então

minimizar

m∑i=1

(fi −M(xi))2



Caso polinomial - estabilidade numérica

No caso polinomial, ou seja, no caso em que M(x) = pn(x) e é usado oconjunto dos polinómios base {1, x, x2, . . . , xn−1, xn} na definição domodelo, temos que

M(x) = pn(x) = c0 + c1x + c2x2 + · · ·+ cn−1x

n−1 + cnxn

e a resolução do problema

minm∑

i=1

(fi − pn(xi))2

resulta num sistema linear mal condicionado.

A introdução dos polinómios ortogonais é suficiente para resolver oproblema de mau condicionamento.



Polinómios ortogonais

Definições:Duas funções g(x) e h(x) dizem-se ortogonais se o seu produtoescalar for nulo, i.e., se 〈g(x), h(x)〉 = 0.A sequência P0(x), P1(x), . . . , Pn(x) forma uma sequência de n + 1polinómios ortogonais se os polinómios Pi(x), i = 0, . . . , n, foremortogonais dois a dois e cada Pi(x) for um polinómio de grau igual a i.



Polinómios ortogonais

Definições:Duas funções g(x) e h(x) dizem-se ortogonais se o seu produtoescalar for nulo, i.e., se 〈g(x), h(x)〉 = 0.A sequência P0(x), P1(x), . . . , Pn(x) forma uma sequência de n + 1polinómios ortogonais se os polinómios Pi(x), i = 0, . . . , n, foremortogonais dois a dois e cada Pi(x) for um polinómio de grau igual a i.



Propriedades dos polinómios ortogonais

Os polinómios ortogonais Pn(x) são também linearmenteindependentes. Assim, qualquer polinómio pn(x), de grau n, pode serexpresso na seguinte forma única

pn(x) = c0P0(x) + c1P1(x) + · · ·+ cnPn(x),

como uma combinação linear de uma sequência de polinómiosortogonais.Os polinómios Pn(x) têm zeros reais e distintos que pertencem a[x1, xm].



Propriedades dos polinómios ortogonais

Os polinómios ortogonais Pn(x) são também linearmenteindependentes. Assim, qualquer polinómio pn(x), de grau n, pode serexpresso na seguinte forma única

pn(x) = c0P0(x) + c1P1(x) + · · ·+ cnPn(x),

como uma combinação linear de uma sequência de polinómiosortogonais.Os polinómios Pn(x) têm zeros reais e distintos que pertencem a[x1, xm].



Determinação dos polinómios ortogonais

Pi+1(x) = Ai(x−Bi)Pi(x)− CiPi−1(x), para i = 0, . . . , n− 1

sendo P−1(x) = 0, P0(x) = 1 e os coeficientes da relação

Ai = 1, para todo o i,

Bi =〈xPi(x), Pi(x)〉〈Pi(x), Pi(x)〉

, para todo o i,

C0 = 0 e Ci =〈Pi(x), Pi(x)〉〈Pi−1(x), Pi−1(x)〉

para i > 0.



Determinação dos coeficientes

c0 =∑m

i=1 fiP0(xi)∑mi=1 P0(xi)2

, c1 =∑m

i=1 fiP1(xi)∑mi=1 P1(xi)2

, . . . ,

cn =∑m

i=1 fiPn(xi)∑mi=1 Pn(xi)2

sendo o modelo

pn(x) = c0P0(x) + c1P1(x) + · · ·+ cnPn(x).



Exemplo

Determinação da recta (polinómio de grau 1) que melhor se ajusta (nosentido dos mínimos quadrados) aos dados da tabela

xi 20 40 50 60 60 70 80fi 5.5 5.6 5.7 5.9 5.85 6.1 7.5

O modelo a determinar será p1(x) = c0P0(x) + c1P1(x). Em primeiro lugardetermina-se os polinómios ortogonais. P0(x) = 1 por definição e

P1(x) = A0(x−B0)P0(x)− C0P−1(x)

com A0 = 1, C0 = 0, P−1(x) = 0 e

B0 =〈xP0(x), P0(x)〉〈P0(x), P0(x)〉

=∑7

i=1 xiP20 (xi)∑7

i=1 P 20 (xi)

=3807

= 54.285714



Cont.

Temos então que p1(x) = c0 + c1(x− 54.285714) e falta determinar oscoeficientes.

c0 =∑7

i=1 fiP0(xi)∑7i=1 P0(xi)2

=42.15

7= 6.021429

c1 =∑7

i=1 fiP1(xi)∑7i=1 P1(xi)2

=5.5(20− 54.285714) · · · 7.5(80− 54.285714)

(20− 54.285714)2 · · · (80− 54.285714)2

=62.5714412371.4286

= 0.026386 (0.026506 CoNum)

ficando p1(x) = 6.021429 + 0.026386(x− 54.285714)



Exemplo - Interpretação gráfica

20 30 40 50 60 70 805

5.5

6

6.5

7

7.5

8

p1(x)

p2(x)

p2(x) = 6.021429 + 0.026506(x− 54.285714)+0.000946 ∗ ((x− 44.991394)(x− 54285714)− 338.775510)



Modelo linear não polinomial - Caso discreto

Pretende-se determinar o modelo, linear nos coeficientes, que melhor seaproxima (no sentido dos mínimos quadrados) à função f(x). Sendo omodelo linear nos coeficientes pode-se escrever como

M(x) = c1φ1(x) + c2φ2(x) + · · ·+ cnφn(x)

ExemploM(x) = c1 + c2e

x

em que φ1(x) = 1 e φ2(x) = ex.



Dedução das equações normais

Pretende-se

minimizarS ≡m∑

i=1

(fi − c1φ1(xi)− · · · − cnφn(xi))2

e derivando S em ordem aos ci, i = 1, . . . , n e igualando a zero, obtém-se

∂S

∂c1=

m∑i=1

(fi − c1φ1(xi)− · · · − cnφn(xi))φ1(xi) = 0

· · ·

∂S

∂cn=

m∑i=1

(fi − c1φ1(xi)− · · · − cnφn(xi))φn(xi) = 0



Sistema das equações normais

Na forma matricial vem0BBBBBBBBBBBB@

mXi=1

φ1(xi)φ1(xi)

mXi=1

φ1(xi)φ2(xi) ...

mXi=1

φ1(xi)φn(xi)

mXi=1

φ2(xi)φ1(xi)

mXi=1

φ2(xi)φ2(xi) ...

mXi=1

φ2(xi)φn(xi)

...... ......mX

i=1

φn(xi)φ1(xi)

mXi=1

φn(xi)φ2(xi) ...

mXi=1

φn(xi)φn(xi)

1CCCCCCCCCCCCA

0BBB@c1

c2

...

cn

1CCCA =

0BBBBBBBBBBBB@

mXi=1

fiφ1(xi)

mXi=1

fiφ2(xi)

...mX

i=1

fiφn(xi)

1CCCCCCCCCCCCA



Exemplo

Consideremos o exemplo já apresentado

xi(velocidade - km/h) 20 40 50 60 60 70 80fi(consumo l/100km) 5.5 5.6 5.7 5.9 5.85 6.1 7.5

pretendendo-se ajustar o modelo M(x) = c1 + c2ex20 à tabela dada. Temos

que φ1(x) = 1 e φ2(x) = ex20 , resolvendo o sistema linear (EGPP)( ∑7

i=1 1∑7

i=1 exi20

∑7i=1 fi∑7

i=1 exi20

∑7i=1 e

xi10

∑7i=1 fie

xi20

)

vem c1 = 5.231132, c2 = 0.036838 e consequentemente

M(x) = 5.231132 + 0.036838ex20 .



Exemplo - Gráfico

10 20 30 40 50 60 70 80 905

5.5

6

6.5

7

7.5

Velocidade (km/h)

Con

sum

o (l/

100k

m)



Resíduo

O resíduo é determinado por

m∑i=1

(fi − c1φ1(xi)− · · · − cnφn(xi))2 .

No caso do polinómio de grau um, p1(x), temos um resíduo de 1.123193.Para o polinómio de grau 2 temos um resíduo de 0.302736 e no modelolinear um resíduo de 0.247494.O melhor modelo é o linear, porque tem menor resíduo.



Exemplo - Gráficos

10 20 30 40 50 60 70 80 905

5.5

6

6.5

7

7.5

Velocidade (km/h)

Con

sum

o (l/

100k

m)

p1(x)

p2(x)

M(x)


Equações diferenciais - Condições iniciais

Conteúdo

1 Introdução

2 Erros












Motivação

As equações diferenciais surgem nas mais diversas áreas da engenharia esão caracterizadas por descrever uma relação entre a variável independente(x) e a dependente (y(x), ou simplesmente y, no casos em que não háambiguidade), envolvendo derivadas.

Forma geralUma equação diferencial é genericamente descrita como

F (x, y(x), y′(x), y′′(x), . . . , y(m−1)(x), y(m)(x)) = 0

ouF (x, y, y′, y′′, . . . , y(m−1), y(m)) = 0



Características

Ordem da equação diferencialA ordem mais elevada da derivada de y que aparece na equação define aordem da equação diferencial.

Exemplo

2y′′ + 3y′ − xy + x− 3 = 0

É uma equação diferencial de ordem 2 em que x é a variável independentee y ≡ y(x) é a variável dependente.



Características

Ordem da equação diferencialA ordem mais elevada da derivada de y que aparece na equação define aordem da equação diferencial.

Exemplo

2y′′ + 3y′ − xy + x− 3 = 0

É uma equação diferencial de ordem 2 em que x é a variável independentee y ≡ y(x) é a variável dependente.



Características

Solução da equação diferencial

A solução de uma equação diferencial é uma função (y(x)) – (Ou famíliade funções).

Exemplo

y′ = sen(x)

Tem como solução y(x) = −cos(x) + c, em que c é uma constante.



Características

Solução da equação diferencial

A solução de uma equação diferencial é uma função (y(x)) – (Ou famíliade funções).

Exemplo

y′ = sen(x)

Tem como solução y(x) = −cos(x) + c, em que c é uma constante.



Características

Equações diferenciais ordináriasUm equação diferencial que possua apenas uma variável independentediz-se ordinária.

Equações diferenciais com derivadas parciaisUma equação diferencial que possua duas ou mais variáveis independentesdiz-se com derivadas parciais.

Exemplo com derivadas parciais

∂2u

∂x2+

∂2u

∂y2= ex+y

É uma equação diferencial de ordem 2 com derivadas parciais (u ≡ u(x, y)).



Características




∂2u

∂x2+

∂2u

∂y2= ex+y




Características




∂2u

∂x2+

∂2u

∂y2= ex+y




Solução única?

Uma equação diferencial só terá solução única se forem especificadascondições auxiliares na função y.

Exemplo

No caso anterior (y′ = sen(x)) se indicarmos que y(0) = 3 temos quey(0) = −cos(0) + c = −1 + c = 3, logo c = 4, ficandoy(x) = −cos(x) + 4 como solução única para a equação diferencial.

Exemplo de segunda ordem

Considere-se a equação diferencial y′′ = x. A solução é (integrando duasvezes) y(x) = (x+a)2

2 + b, sendo a e b duas constantes.

Regra geralUma equação diferencial de ordem m, para ter solução única, necessita dem condições auxiliares.



Solução única?


Exemplo








Solução única?


Exemplo








Resolução

Na prática a determinação de uma solução analítica pode não ser possível(é mesmo o mais provável quando se trata de modelos de casos reais).

A solução numérica consiste numa sequência de valores

y(x0), y(x1), . . . , y(xn−1), y(xn)

que são aproximações à função y(x) nos pontos x0, x1, . . . , xn−1 e xn.

Os pontos x0, x1, . . . , xn−1 e xn são, usualmente, igualmente espaçadosde espaçamento constante e igual a h.

Fala-se pois de uma discretização da equação diferencial porque sedetermina a solução para um conjunto discreto de valores da função.



Resolução



y(x0), y(x1), . . . , y(xn−1), y(xn)






Resolução



y(x0), y(x1), . . . , y(xn−1), y(xn)






Resolução



y(x0), y(x1), . . . , y(xn−1), y(xn)






Equações diferenciais com condições auxiliares

Condições iniciaisQuando as m condições auxiliares são definidas no início do intervalo depontos, para os quais se pretende determinar o valor da função, dizemosque estamos na presença de um problema com condições iniciais.

Condições fronteiraQuando as m condições auxiliares são definidas no início e fim do intervalode pontos, para os quais se pretende determinar o valor da função, dizemosque estamos na presença de um problema com condições fronteira.



Equações diferenciais com condições auxiliares

Condições iniciaisQuando as m condições auxiliares são definidas no início do intervalo depontos, para os quais se pretende determinar o valor da função, dizemosque estamos na presença de um problema com condições iniciais.

Condições fronteiraQuando as m condições auxiliares são definidas no início e fim do intervalode pontos, para os quais se pretende determinar o valor da função, dizemosque estamos na presença de um problema com condições fronteira.



Exemplo de condições auxiliares

Exemplos

Considere-se uma equação diferencial de ordem 2 (y′′ = . . . ) em que sepretende determinar a solução numérica para os pontos 0.1, 0.2, 0.3, 0.4 e0.5.Estaremos na presença de uma equação diferencial com condições iniciaisse, por exemplo, tivermos y(0.1) = 1 e y′(0.1) = 1 como condiçõesauxiliares.Estaremos na presença de uma equação diferencial com condições fronteirase, por exemplo, tivermos y(0.1) = 1 e y′(0.5) = 1 como condiçõesauxiliares.

Condições auxiliaresOs diferentes tipos de condições auxiliares afectam de forma drástica ométodo a aplicar.



Exemplo de condições auxiliares

Exemplos

Considere-se uma equação diferencial de ordem 2 (y′′ = . . . ) em que sepretende determinar a solução numérica para os pontos 0.1, 0.2, 0.3, 0.4 e0.5.Estaremos na presença de uma equação diferencial com condições iniciaisse, por exemplo, tivermos y(0.1) = 1 e y′(0.1) = 1 como condiçõesauxiliares.Estaremos na presença de uma equação diferencial com condições fronteirase, por exemplo, tivermos y(0.1) = 1 e y′(0.5) = 1 como condiçõesauxiliares.

Condições auxiliaresOs diferentes tipos de condições auxiliares afectam de forma drástica ométodo a aplicar.



Equações diferenciais ordinárias de 1a ordem com condiçõesiniciais

Forma geral

y′ = f(x, y) y(x0) = y0

Método de passo único. Fórmula de Euler e de Runge-Kutta de 1a ordem

yi+1 = yi + hf(xi, yi), i = 0, 1, . . . , n− 1

Designámos de etapa a determinação do valor de yi+1 à custa dos valoresde yi e xi.




Forma geral

y′ = f(x, y) y(x0) = y0


yi+1 = yi + hf(xi, yi), i = 0, 1, . . . , n− 1





Forma geral

y′ = f(x, y) y(x0) = y0


yi+1 = yi + hf(xi, yi), i = 0, 1, . . . , n− 1




Exemplo

Exemplo

Considere-se a equação diferencial y′ + y = sen(x)− cos(x) com acondição inicial de y(0) = −1. Pretende-se determinar os valores de y(x)para 0.15, 0.30, 0.45, 0.60 através da fórmula de Runge-Kutta de 1a

ordem.

Dados do problema

Da equação diferencial temos que f(x, y) = sen(x)− cos(x)− y.h = 0.15 e são necessárias 4 etapas.

i 0 1 2 3 4xi 0.00 0.15 0.30 0.45 0.60yi −1.0000 ? ? ? ?



Exemplo

Exemplo

Considere-se a equação diferencial y′ + y = sen(x)− cos(x) com acondição inicial de y(0) = −1. Pretende-se determinar os valores de y(x)para 0.15, 0.30, 0.45, 0.60 através da fórmula de Runge-Kutta de 1a

ordem.

Dados do problema

Da equação diferencial temos que f(x, y) = sen(x)− cos(x)− y.h = 0.15 e são necessárias 4 etapas.

i 0 1 2 3 4xi 0.00 0.15 0.30 0.45 0.60yi −1.0000 ? ? ? ?



Etapas

1a etapa

y1 = y0 + hf(x0, y0) = y0 + h(sen(x0)− cos(x0)− y0) =− 1 + 0.15(sen(0)− cos(0) + 1) = −1

2a etapa

y2 = −1 + 0.15(sen(0.15)− cos(0.15) + 1) = −0.9759

3a etapa

y3 = −0.9759 + 0.15(sen(0.30)− cos(0.30) + 0.9759) = −0.9285

4a etapa

y4 = −0.9285 + 0.15(sen(0.45)− cos(0.45) + 0.9285) = −0.8590



Gráfico

0 0.1 0.2 0.3 0.4 0.5 0.6−1

−0.98

−0.96

−0.94

−0.92

−0.9

−0.88

−0.86

−0.84

Condição inicial Erro de truncatura



Erro de truncatura

Erro de truncaturaO erro de truncatura é dado por

eT =12h2y′′(ξi) para ξi ∈ [xi, xi+1].

Resultante da expansão de y em série de Taylor

y(x + h) = y(x) + h ∗ y′(x) +12h2y′′(ξ).



Exemplo

Cálculo do erroO erro de truncatura pode ser calculado implicitamente através do uso daprópria equação diferencial.

y′′ = (y′)′ = (sen(x)− cos(x)− y)′ = cos(x) + sen(x)− y′ =cos(x) + sen(x)− sen(x) + cos(x) + y = 2cos(x) + y

Um majorante para y′′ pode ser 2 ∗ cos(0.0) + 1 = 3. Resultando nummajorante para o valor absoluto do erro de truncatura de

12× 0.152 ∗ 3 = 0.03375



Runge-Kutta de segunda ordem

Fórmula

yi+1 = yi +12(p + q), i = 1, 2, . . .

comp = hf(xi, yi) e q = hf(xi+1, yi + p)

Exemplo y′ = sen(x)− cos(x)− y, y(0) = −1

p =hf(x0, y0) = 0.15(sen(0)− cos(0) + 1) = 0q =hf(x1, y0 + p) = hf(0.15,−1 + 0) =

0.15(sen(0.15)− cos(0.15) + 1) = 0.0241

resultado emy1 = y0 +

12(p + q) = −1 +

12(0 + 0.0241) = −0.9880



Runge-Kutta de segunda ordem

Fórmula

yi+1 = yi +12(p + q), i = 1, 2, . . .

comp = hf(xi, yi) e q = hf(xi+1, yi + p)

Exemplo y′ = sen(x)− cos(x)− y, y(0) = −1

p =hf(x0, y0) = 0.15(sen(0)− cos(0) + 1) = 0q =hf(x1, y0 + p) = hf(0.15,−1 + 0) =

0.15(sen(0.15)− cos(0.15) + 1) = 0.0241

resultado emy1 = y0 +

12(p + q) = −1 +

12(0 + 0.0241) = −0.9880



Gráfico

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7−1

−0.98

−0.96

−0.94

−0.92

−0.9

−0.88

−0.86

−0.84

−0.82

x

y



Sistemas de equações diferenciais de primeira ordem –Runge-Kutta de 2a ordem

Forma geraly′1(x) = f1(x, y1(x), . . . , yk(x))y′2(x) = f2(x, y1(x), . . . , yk(x))

. . .y′k(x) = fk(x, y1(x), . . . , yk(x))

y1(x0) = y10

y2(x0) = y20

. . .yk(x0) = yk0

Existem k equações diferenciais de 1a ordem que relacionam k funções y.

Exemplo

{y′1(x) = f1(x, y1(x), y2(x)) = y2(x)− x2

y′2(x) = f2(x, y1(x), y2(x)) = 2x− y1(x)

{y1(0.5) = 0.4794y2(0.5) = 1.1276

Cuja solução analítica é y1(x) = sen(x) y2(x) = cos(x) + x2



Sistemas de equações diferenciais de primeira ordem –Runge-Kutta de 2a ordem

Forma geraly′1(x) = f1(x, y1(x), . . . , yk(x))y′2(x) = f2(x, y1(x), . . . , yk(x))

. . .y′k(x) = fk(x, y1(x), . . . , yk(x))

y1(x0) = y10

y2(x0) = y20

. . .yk(x0) = yk0

Existem k equações diferenciais de 1a ordem que relacionam k funções y.

Exemplo

{y′1(x) = f1(x, y1(x), y2(x)) = y2(x)− x2

y′2(x) = f2(x, y1(x), y2(x)) = 2x− y1(x)

{y1(0.5) = 0.4794y2(0.5) = 1.1276

Cuja solução analítica é y1(x) = sen(x) y2(x) = cos(x) + x2



Equações

Equações do método de Runge-Kutta de 2a ordem, i = 0, 1, . . .

p1 = hf1(xi, y1i, . . . , yki). . .

pk = hfk(xi, y1i, . . . , yki)

q1 = hf1(xi+1, y1i + p1, . . . , yki + pk). . .

qk = hfk(xi+1, y1i + p1, . . . , yki + pk)

y1,i+1 = y1i +12(p1 + q1)

. . .

yk,i+1 = yki +12(pk + qk)



Exemplo – Uma etapa

Notação simplificada{y′1 = y2 − x2

y′2 = 2x− y1

{y10 = 0.4794y20 = 1.1276

x0 = 0.5h = 0.1

p1 =hf1(x0, y10, y20) = 0.1(1.1276− 0.52) = 0.0878p2 =hf2(x0, y10, y20) = 0.1(2× 0.5− 0.4794) = 0.0521q1 =hf1(x1, y10 + p1, y20 + p2)

=0.1f1(0.6, 0.4794 + 0.0878, 1.1276 + 0.0521)=0.1f1(0.6, 0.5672, 1.1797) = 0.0820

q2 =hf2(x1, y10 + p1, y20 + p2) = 0.0633

y11 =y10 +12(p1 + q1) = 0.4794 +

12(0.0878 + 0.0820) = 0.5643

y21 =y20 +12(p2 + q2) = 1.1276 +

12(0.0521 + 0.0633) = 1.1853



Exemplo – Uma etapa

Notação simplificada{y′1 = y2 − x2

y′2 = 2x− y1

{y10 = 0.4794y20 = 1.1276

x0 = 0.5h = 0.1

p1 =hf1(x0, y10, y20) = 0.1(1.1276− 0.52) = 0.0878p2 =hf2(x0, y10, y20) = 0.1(2× 0.5− 0.4794) = 0.0521q1 =hf1(x1, y10 + p1, y20 + p2)

=0.1f1(0.6, 0.4794 + 0.0878, 1.1276 + 0.0521)=0.1f1(0.6, 0.5672, 1.1797) = 0.0820

q2 =hf2(x1, y10 + p1, y20 + p2) = 0.0633

y11 =y10 +12(p1 + q1) = 0.4794 +

12(0.0878 + 0.0820) = 0.5643

y21 =y20 +12(p2 + q2) = 1.1276 +

12(0.0521 + 0.0633) = 1.1853



Gráfico

0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 10.4

0.6

0.8

1

1.2

1.4

1.6

1.8

x

y 1(x),

y 2(x)

y1=sen(x)

y2=cos(x)+x2

y20

=y(0.5)

y10

y23

y11

=y(0.6)

y25

=y(1.0)

x0=0.5

h=0.1



Gráfico – Zoom para verificação de erro

0.6975 0.698 0.6985 0.699 0.6995 0.7 0.7005 0.701 0.7015 0.702

0.638

0.64

0.642

0.644

0.646

0.648

0.65

x

y 1(x)

y13



Equação diferencial de ordem superior com condições iniciais

Forma geral

y(n)(x) = f(x, y(x), y′(x), . . . , y(n−1)(x))

y(x0) = c0, y′(x0) = c1, . . . , y(n−1)(x0) = cn−1

Exemplo

y′′′(x)+y′′(x)+y′(x)+y(x) = 0 em que y(0) = −1, y′(0) = 0, y′′(0) = 1

i.e.y′′′ = −y′′ − y′ − y cuja solução analítica é y = −cos(x)



Equação diferencial de ordem superior com condições iniciais

Forma geral

y(n)(x) = f(x, y(x), y′(x), . . . , y(n−1)(x))

y(x0) = c0, y′(x0) = c1, . . . , y(n−1)(x0) = cn−1

Exemplo

y′′′(x)+y′′(x)+y′(x)+y(x) = 0 em que y(0) = −1, y′(0) = 0, y′′(0) = 1

i.e.y′′′ = −y′′ − y′ − y cuja solução analítica é y = −cos(x)



Abordagem

Transformação em sistema de equações de primeira ordemIntrodução de n variáveis auxiliares yi, i = 1, . . . , n e escrita do sistema

y′1 = y2

y′2 = y3

y′3 = y4

. . .y′n−1 = yn

yn = f(x, y1, y2, . . . , yn−1, yn)

e as respectivas condições iniciais

y10 = c0, y20 = c1, . . . , yn0 = cn−1



Exemplo

Exemplo

y′′′ = −y′′ − y′ − y, y(0) = −1, y′(0) = 0, y′′(0) = 1

Temos que n = 3, introduzindo as 3 variáveis adicionais y1 = y, y2 = y′ ey3 = y′′ vem que:

y′1 = y2

y′2 = y3

y′3 = −y3 − y2 − y1

Com as seguintes condições iniciais

y1(0) = y10 = −1, y2(0) = y20 = 0, y3(0) = y30 = 1



Uma etapa com h = 0.1

A resolução do sistema pelo método de Runge-Kutta de segunda ordemresulta em:

p1 =hf1(x0, y10, y20, y30) = 0.1f1(0,−1, 0, 1) = 0.1× 0 = 0p2 =hf2(x0, y10, y20, y30) = 0.1× 1 = 0.1p3 =hf3(x0, y10, y20, y30) = 0.1× (−1− 0 + 1) = 0q1 =hf1(x1, y10 + p1, y20 + p2, y30 + p3)

=0.1f1(0.1,−1, 0.1, 1) = 0.1× 0.1 = 0.01q2 =0.1f2(0.1,−1, 0.1, 1) = 0.1× 1 = 0.1q3 =0.1f3(0.1,−1, 0.1, 1) = 0.1(−1− 0.1 + 1) = −0.01

y11 =y10 +12(p1 + q1) = −1 +

12(0 + 0.01) = −0.9950

y21 =y20 +12(p2 + q2) = 0 +

12(0.1 + 0.1) = 0.1

y31 =y30 +12(p2 + q3) = 1 +

12(0− 0.01) = 0.9950



Gráfico

0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5−1

−0.8

−0.6

−0.4

−0.2

0

0.2

0.4

0.6

0.8

1

x

y(x)

, y’(x

), y

’’(x)

y2(x)=y’(x)=sen(x)

y3(x)=y’’(x)=cos(x)

y1(x)=y(x)=−cos(x)



Teste de conhecimentos

Questão

Seríamos capazes de obter y(4)(0.1) com os dados obtidos até agora?

RespostaClaro que sim...

y(4) = (y′′′)′ = (−y′′−y′−y)′ = −y′′′−y′′−y′ = −(−y′′−y′−y)−y′′−y′ = y

y(4) = y, logo y(4)(0.1) = y(0.1) = y11 = −0.9950



Teste de conhecimentos

Questão

Seríamos capazes de obter y(4)(0.1) com os dados obtidos até agora?

RespostaClaro que sim...

y(4) = (y′′′)′ = (−y′′−y′−y)′ = −y′′′−y′′−y′ = −(−y′′−y′−y)−y′′−y′ = y

y(4) = y, logo y(4)(0.1) = y(0.1) = y11 = −0.9950



Exemplo mais complexo

Exemplo com duas equações (uma de ordem superior){u′′ + u′v′ = sen(x)v′ + v + u = cos(x)

→{

u′′ = sen(x)− u′v′

v′ = cos(x)− v − u

com x ∈ [0, 1], u(0) = u′(0) = 0 e v(0) = 1.

ResoluçãoPara a primeira equação serão necessárias duas variáveis auxiliares e comotal u→ y1, v → y3. Fica então

y′1 = y2

y′2 = sen(x)− y2y′3

y′3 = cos(x)− y3 − y1

→

y′1 = y2

y′2 = sen(x)− y2(cos(x)− y3 − y1)y′3 = cos(x)− y3 − y1

com x ∈ [0, 1], y1(0) = y2(0) = 0 e y3(0) = 1.Aplicação do método de Runge-Kutta para o sistema de três equações.



Exemplo mais complexo

Exemplo com duas equações (uma de ordem superior){u′′ + u′v′ = sen(x)v′ + v + u = cos(x)

→{

u′′ = sen(x)− u′v′

v′ = cos(x)− v − u

com x ∈ [0, 1], u(0) = u′(0) = 0 e v(0) = 1.

ResoluçãoPara a primeira equação serão necessárias duas variáveis auxiliares e comotal u→ y1, v → y3. Fica então

y′1 = y2

y′2 = sen(x)− y2y′3

y′3 = cos(x)− y3 − y1

→

y′1 = y2

y′2 = sen(x)− y2(cos(x)− y3 − y1)y′3 = cos(x)− y3 − y1

com x ∈ [0, 1], y1(0) = y2(0) = 0 e y3(0) = 1.Aplicação do método de Runge-Kutta para o sistema de três equações.


Equações diferenciais - Condições Fronteira

Conteúdo

1 Introdução

2 Erros












Equações diferenciais com condições fronteira


F (x, y(x), y′(x), y′′(x), . . . , y(m−1)(x), y(m)(x)) = 0

ouF (x, y, y′, y′′, . . . , y(m−1), y(m)) = 0

Condições fronteiraTrata-se de determinar a solução num conjunto de pontos x0, x1, ..., xn

em que são conhecidas m condições auxiliares envolvendo y, y′, . . . , ym

em x0 e xn.



Equações diferenciais com condições fronteira


F (x, y(x), y′(x), y′′(x), . . . , y(m−1)(x), y(m)(x)) = 0

ouF (x, y, y′, y′′, . . . , y(m−1), y(m)) = 0

Condições fronteiraTrata-se de determinar a solução num conjunto de pontos x0, x1, ..., xn

em que são conhecidas m condições auxiliares envolvendo y, y′, . . . , ym

em x0 e xn.



Exemplo

Exemplo de equação de segunda ordem

y′′ + y′ + y = cos(x)

com as seguintes condições auxiliares (fronteira)

y′(0) + y′(1) = 1.5403 Qualquer relação linear de y, y′, y′′

y′(0) + y(1) = 1.8415 nos pontos fronteira do intervalo

Pretende-se determinar a solução no intervalo [0, 1] com h = 0.25.Solução analítica y = sen(x).



Técnica de resolução

ResoluçãoUma das possíveis técnicas de resolução consiste no uso de diferençasfinitas para aproximar as derivadas envolvidas na equações.Todos os valores de y são calculados em simultâneo através da resoluçãode um sistema linear.



Diferenças finitas

Exemplo de dedução das diferenças finitas

y(xi+1) = y(xi + h) = y(xi) + hy′(xi) + 12h2y′′(ξ) ξ ∈ [xi, xi+1]

- y(xi−1) = y(xi − h) = y(xi)− hy′(xi) + 12h2y′′(η) η ∈ [xi−1, xi]

y(xi+1)− y(xi−1) = 2hy′(xi) + 12h2y′′(ξ)− 1

2h2y′′(η)

Diferença dividida central de primeira ordem

y′(xi) ≈(y(xi+1)− y(xi−1))

2h=

yi+1 − yi−1

2h

Diferença dividida central porque relaciona a derivada em xi com a funçãoem xi+1 e xi−1.



Diferenças finitas




y(xi+1)− y(xi−1) = 2hy′(xi) + 12h2y′′(ξ)− 1

2h2y′′(η)


y′(xi) ≈(y(xi+1)− y(xi−1))

2h=

yi+1 − yi−1

2h




Diferenças finitas




y(xi+1)− y(xi−1) = 2hy′(xi) + 12h2y′′(ξ)− 1

2h2y′′(η)


y′(xi) ≈(y(xi+1)− y(xi−1))

2h=

yi+1 − yi−1

2h




Diferenças finitas centrais e ascendentes

Centrais

y′i ≈yi+1 − yi−1

2h

y′′i ≈yi+1 − 2yi + yi−1

h2

y′′′i ≈yi+2 − 2yi+1 + 2yi−1 − yi−2

2h3

Ascendentes

y′i ≈yi − yi−1

h

y′′i ≈yi − 2yi−1 + yi−2

h2

y′′′i ≈yi − 3yi−1 + 3yi−2 − yi−3

h3



Diferenças finitas centrais e ascendentes

Centrais

y′i ≈yi+1 − yi−1

2h

y′′i ≈yi+1 − 2yi + yi−1

h2

y′′′i ≈yi+2 − 2yi+1 + 2yi−1 − yi−2

2h3

Ascendentes

y′i ≈yi − yi−1

h

y′′i ≈yi − 2yi−1 + yi−2

h2

y′′′i ≈yi − 3yi−1 + 3yi−2 − yi−3

h3



Diferenças finitas descendentes

Descendentes

y′i ≈yi+1 − yi

h

y′′i ≈yi+2 − 2yi+1 + yi

h2

y′′′i ≈yi+3 − 3yi+2 + 3yi+1 − yi

h3



Resolução exemplo

Discretização (intervalo [0, 1] com h = 0.25)

i 0 1 2 3 4xi 0 0.25 0.5 0.75 1yi ? ? ? ? ?

Na equação diferencial usam-se as diferenças finitas centrais. Temos poispara cada valor de i (i = 1, 2, 3) na equação diferencial:

Substituição

y′′i + y′i + yi = cos(xi), i = 1, 2, 3

i.e.

yi+1 − 2yi + yi−1

h2+

yi+1 − yi−1

2h+ yi = cos(xi), i = 1, 2, 3



Resolução exemplo


i 0 1 2 3 4xi 0 0.25 0.5 0.75 1yi ? ? ? ? ?


Substituição

y′′i + y′i + yi = cos(xi), i = 1, 2, 3

i.e.

yi+1 − 2yi + yi−1

h2+

yi+1 − yi−1

2h+ yi = cos(xi), i = 1, 2, 3



Resolução exemplo


i 0 1 2 3 4xi 0 0.25 0.5 0.75 1yi ? ? ? ? ?


Substituição

y′′i + y′i + yi = cos(xi), i = 1, 2, 3

i.e.

yi+1 − 2yi + yi−1

h2+

yi+1 − yi−1

2h+ yi = cos(xi), i = 1, 2, 3



Resolução – Continuação

NotaAs incógnitas são yi, i = 0, 1, 2, 3, 4

Antes de procedermos à substituição de i por 1, 2, 3 vamos simplificar emordem às incógnitas, substituindo h por 0.25 (este processo pode serdiferente, dependendo da complexidade da expressão).

16(yi+1 − 2yi + yi−1) + 2(yi+1 − yi−1) + yi = cos(xi)18yi+1 − 31yi + 14yi−1 = cos(xi)




NotaAs incógnitas são yi, i = 0, 1, 2, 3, 4

Antes de procedermos à substituição de i por 1, 2, 3 vamos simplificar emordem às incógnitas, substituindo h por 0.25 (este processo pode serdiferente, dependendo da complexidade da expressão).

16(yi+1 − 2yi + yi−1) + 2(yi+1 − yi−1) + yi = cos(xi)18yi+1 − 31yi + 14yi−1 = cos(xi)




Discretização

(i = 1) 18y2 − 31y1 + 14y0 = cos(x1) = 0.9689(i = 2) 18y3 − 31y2 + 14y1 = cos(x2) = 0.8776(i = 3) 18y4 − 31y3 + 14y2 = cos(x3) = 0.7317

Três equações lineares a cinco incógnitas. São necessárias as duascondições adicionais.

Condições auxiliares (fronteira)

y′0 + y′4 = 1.5403y′0 + y4 = 1.8415




Discretização

(i = 1) 18y2 − 31y1 + 14y0 = cos(x1) = 0.9689(i = 2) 18y3 − 31y2 + 14y1 = cos(x2) = 0.8776(i = 3) 18y4 − 31y3 + 14y2 = cos(x3) = 0.7317

Três equações lineares a cinco incógnitas. São necessárias as duascondições adicionais.

Condições auxiliares (fronteira)

y′0 + y′4 = 1.5403y′0 + y4 = 1.8415



Resolução – Condições fronteira

Usam-se as diferenças finitas para substituir nas equações auxiliares.

y′(0) + y′(1) =y1 − y0

h+

y4 − y3

h= 1.5403

(Descendente + Ascendente)4y1 − 4y0 + 4y4 − 4y3 = 1.5403

y′(0) + y(1) =y1 − y0

h+ y4 = 1.8415

(Descendente)4y1 − 4y0 + y4 = 1.8415



Resolução – Sistema linear

Considerando as cinco equações resulta no seguinte sistema linear

Sistema linear −4 4 0 −4 4 1.540314 −31 18 0 0 0.96890 14 −31 18 0 0.87760 0 14 −31 18 0.7317−4 4 0 0 1 1.8415

cuja solução (por EGPP) é

Solução y0 = y(0.00) = 0.0936y1 = y(0.25) = 0.3317y2 = y(0.50) = 0.5524y3 = y(0.75) = 0.7420y4 = y(1.00) = 0.8890



Resolução – Sistema linear

Considerando as cinco equações resulta no seguinte sistema linear

Sistema linear −4 4 0 −4 4 1.540314 −31 18 0 0 0.96890 14 −31 18 0 0.87760 0 14 −31 18 0.7317−4 4 0 0 1 1.8415

cuja solução (por EGPP) é

Solução y0 = y(0.00) = 0.0936y1 = y(0.25) = 0.3317y2 = y(0.50) = 0.5524y3 = y(0.75) = 0.7420y4 = y(1.00) = 0.8890



Gráfico

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x

y(x)

y(x)

y2



Gráfico – h = 0.01 – Sistema linear de 101× 101

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 10

0.1

0.2

0.3

0.4

0.5

0.6

0.7

0.8

0.9

x

y(x)

h=0.25

h=0.01



Árvore de decisão

Equação diferencial

Fronteira

Discretização daequação diferencial.Discretização dascondições fron-teira. Resolução desistema linear.

Iniciais

1a ordem Sistema deequações de1a ordem

Ordem supe-rior

Euler ouRunge-Kutta

Runge-Kutta Transformaçãoem sistemade equaçõeslineares.Runge-Kutta.



Mais um exemplo

EnunciadoA equação que representa a perda de calor, T , por uma aleta (secçãoarticulada perto da borda anterior de cada asa de um avião) é a seguinte:

d2T

dx2− aT = −aTambiente

com as seguintes condições:

T (0) = Tparede edT (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. Paraa = 20m−2oC−1, L = 0.3m, Tparede = 200oC e Tambiente = 20oC, estimea temperatura em quatro pontos igualmente distanciados da aleta(considerando o zero como ponto a estimar).



Resolução

Substituindo T por y e as constantes na equação diferencial e nascondições auxiliares obtemos:

y′′ − 20y = −400 com y0 = y(0) = 200 e y′3 = y′(0.3) = 0

i 0 1 2 3xi 0 0.1 0.2 0.3yi 200 ? ? ?

Apesar de conhecermos o valor de y(0) trata-se de uma problema comcondições fronteira, uma vez que a segunda condição auxiliar relaciona aderivada em 0.3.



Discretização

A equação discretizada fica:

yi+1 − 2yi + yi−1

h2−20yi = −400 i.e. 100yi+1−220yi+100yi−1 = −400

Para i = 1, 2 fica:{100y2 − 220y1 = −400− 100× 200 = −20400

100y3 − 220y2 + 100y1 = −400 = −400

A segunda condição resulta em

y′3 =y3 − y2

h= 10y3 − 10y2 = 0



Discretização


yi+1 − 2yi + yi−1

h2−20yi = −400 i.e. 100yi+1−220yi+100yi−1 = −400

Para i = 1, 2 fica:{100y2 − 220y1 = −400− 100× 200 = −20400

100y3 − 220y2 + 100y1 = −400 = −400


y′3 =y3 − y2

h= 10y3 − 10y2 = 0



Discretização


yi+1 − 2yi + yi−1

h2−20yi = −400 i.e. 100yi+1−220yi+100yi−1 = −400

Para i = 1, 2 fica:{100y2 − 220y1 = −400− 100× 200 = −20400

100y3 − 220y2 + 100y1 = −400 = −400


y′3 =y3 − y2

h= 10y3 − 10y2 = 0



Sistema linear

O sistema linear a resolver é então−220 100 0 −20400

100 −220 100 −4000 −10 10 0

y1 y2 y3

cuja solução é

y1 = 151.707317y2 = 129.756098y3 = 129.756098

e y0 = 200 já conhecido.


Integração numérica

Conteúdo

1 Introdução

2 Erros













Pretende-se calcular uma aproximação ao integral definido

I =∫ b

af(x)dx

com a e b constantes.A técnica a utilizar consiste na aproximação da função f(x) por umpolinómio interpolador, pn(x), e posteriormente aproxima-se I por∫ ba pn(x)dx.



Motivação para a integração numérica

As situações mais frequentes onde se torna necessário calcular umaaproximação ao integral definido são:

a função primitiva não pode vir expressa em termos de funçõeselementares;a expressão da primitiva é muito complexa;a função integranda é conhecida apenas num conjunto discreto depontos.













Fórmulas com intervalos de amplitudes constantes

Fórmulas Newton-Cotes (simples)As fórmulas simples são aplicadas a um intervalo que contém o númeromínimo de pontos para a interpolação entre a e b.Regra do Trapézio (n = 1)∫ b

af(x)dx ≈ b− a

2[f(a) + f(b)]



Cont.

Regra de Simpson (n = 2)∫ b

af(x)dx ≈ b− a

6

[f(a) + 4f

(a + b

2

)+ f(b)

]

Regra dos 38 (n = 3, regra de Newton dos três oitavos)∫ b

af(x)dx ≈ b− a

8

[f(a) + 3f

(2a + b

3

)+ 3f

(a + 2b

3

)+ f(b)

]



Interpretação geométrica - Trapézio

0 0.5 1 1.5 2 2.5 3 3.50

0.2

0.4

0.6

0.8

1

1.2

1.4

x

sin(

x)

p1(x)

sin(x)

I

p1(x) =sen(1) + sen(1)−sen(1.5)

1−1.5 (x− 1)



Erro de truncatura

As fórmulas do erro de truncatura são deduzidas a partir da fórmula para oerro da interpolação polinomial.Regra do Trapézio

ETT = −(b− a)3

12f ′′(ξ), ξ ∈ [a, b]

Uma vez que ξ é um valor desconhecido calcula-se um majorante do erroabsoluto, através da expressão∣∣∣∣−(b− a)3

12M

∣∣∣∣ , em que M ≥ maxξ∈[a,b]

f ′′(ξ)



Cont.

Regra de Simpson

ETS = −(b− a)5

2880f (iv)(ξ), ξ ∈ [a, b]

Regra dos 38

ET 38

= −(b− a)5

6480f (iv)(ξ), ξ ∈ [a, b]



Fórmulas compostas - caso do Trapézio

Suponhamos que se pretende aplicar a regra do Trapézio a um conjunto depontos, igualmente espaçados, x0, x1, x2, . . . , xn, para determinar∫ xn

x0f(x)dx (dados os valores de f(x0) = f0, . . . , f(xn) = fn).

Como∫ xn

x0

f(x)dx =∫ x1

x0

f(x)dx +∫ x2

x1

f(x)dx + · · ·+∫ xn

xn−1f(x)dx

e ∫ x1

x0

f(x)dx ≈ x1 − x0

2[f0 + f1] ,∫ xi+1

xi

f(x)dx ≈ xi+1 − xi

2[fi + fi+1]



Cont.

Temos então que∫ xn

x0

f(x)dx ≈ x1 − x0

2[f0 + f1] +

x2 − x1

2[f1 + f2] + . . .

+xn − xn−1

2[fn−1 + fn]

Usando o facto de que h = x1 − x0 = x2 − x1 = · · · = xn − xn−1, ecolocando em evidencia o termo h

2 vem que∫ xn

x0

f(x)dx ≈ h

2[f0 + 2f1 + 2f2 + · · ·+ 2fn−1 + fn]



Cont.

As fórmulas compostas de Simpson e 38 são deduzidas de forma semelhante

à do Trapézio, mas no caso da regra de Simpson consideram-se múltiplosde 2 intervalos e no caso dos 3

8 múltiplos de 3 intervalos.

Composta de Simpson∫ xn

x0

f(x)dx ≈ h

3[f0 + 4f1 + 2f2 + 4f3 + · · ·+ 2fn−2 + 4fn−1 + fn]

Composta dos 38∫ xn

x0

f(x)dx ≈ 3h

8[f0 + 3f1 + 3f2 + 2f3 + · · ·+

+ · · ·+ 2fn−3 + 3fn−2 + 3fn−1 + fn]



Erros de truncatura

As fórmulas dos erros de truncatura são deduzidas a partir da fórmula doerro de truncatura simples.Trapézio

ETT (h) =−n−1∑i=0

(xi+1 − xi)3

12f ′′(ξi) = −h2

12(b− a)

n

n−1∑i=0

f ′′(ξi)

=− h2

12(b− a)f ′′(η), η ∈ [a, b].



Cont.

Simpson

ETS(h) = − h4

180(b− a)f (iv)(η), η ∈ [a, b]

38

ET 38(h) = −h4

80(b− a)f (iv)(η), η ∈ [a, b]



Intervalos de amplitudes variadas

Os intervalos de amplitudes variadas devem ser divididos em intervalos deamplitudes constantes e posteriormente aplicar-se as fórmulasNewton-Cotes.

Da análise das fórmulas do erro de truncatura conclui-se a aplicação dasseguintes regras (para valores de h ≤ 1):

quando o número de intervalos é par aplica-se a regra de Simpson;quando o número de intervalos é múltiplo de 3 aplica-se a regra dos 3

8 ;aplica-se a regra do trapézio nos restantes casos.

















Exemplo [Edite Fernandes, 1998, adaptado]

Um carro inicia um percurso, num dia de Inverno, e um aparelho mede oconsumo de gasolina no instante em que percorreu x km. Os resultadosregistados são:

x km 0.00 1.25 2.50 3.75 5.00 6.50 8.00 9.50 10.00f(x) l/hm 0.260 0.208 0.172 0.145 0.126 0.113 0.104 0.097 0.092

em que f(x) designa o consumo no fim do percurso, naquele instante, eml/km. Calcule o consumo total de gasolina no fim do percurso de 10km.∫ 10

0f(x)dx ≈1.25

3[0.260 + 4× 0.208 + 2× 0.172 + 4× 0.145 + 0.126]

+3× 1.5

8[0.126 + 3× 0.113 + 3× 0.104 + 0.097]

+0.52

[0.097 + 0.092]



Erro absoluto de truncatura

É a soma dos majorantes das fórmulas aplicadas. Como a função f(x) édesconhecida usamos uma tabela das diferenças divididas para as estimar.

ET =∣∣∣∣−1.254

180(5.00− 0.00)0.000102× 4!

∣∣∣∣+∣∣∣∣−1.54

80(9.50− 5.00)0.000016× 4!

∣∣∣∣+∣∣∣∣−0.52

12(10− 9.5)(−0.002667)× 2!

∣∣∣∣



Cont.

x f(x)0.00 0.260

-0.0416001.25 0.208 0.005120

-0.028800 -0.0005972.50 0.172 0.002880 0.000102

-0.021600 -0.0000853.75 0.145 0.002560 0.000007

-0.015200 -0.0000465.00 0.126 0.002376 -0.000055

-0.008667 -0.0003506.50 0.113 0.000889 0.000044

-0.006000 -0.0000998.00 0.104 0.000444 0.000016

-0.004667 -0.0000009.50 0.097 -0.002667

-0.01000010.00 0.092



Exercício

EnunciadoConsidere o seguinte integral

I =∫ 1

0xexdx.

Calcule uma aproximação ao integral I obtida a partir da fórmula compostade Simpson, de tal forma que o erro (de truncatura) cometido, em valorabsoluto, não exceda 0.005.

Temos então que

f(x) = xex, f ′(x) = (x + 1)ex, f ′′(x) = (x + 2)ex,

f ′′′(x) = (x + 3)ex, f (iv)(x) = (x + 4)ex



Exercício

EnunciadoConsidere o seguinte integral

I =∫ 1

0xexdx.

Calcule uma aproximação ao integral I obtida a partir da fórmula compostade Simpson, de tal forma que o erro (de truncatura) cometido, em valorabsoluto, não exceda 0.005.

Temos então que

f(x) = xex, f ′(x) = (x + 1)ex, f ′′(x) = (x + 2)ex,

f ′′′(x) = (x + 3)ex, f (iv)(x) = (x + 4)ex



Exercício

Majorante

O majorante de f (iv) no intervalo [0, 1] é então

M = majη∈[0,1]f(iv)(η) = (1 + 4)e1 = 5e

Erro de truncaturaDa fórmula composta de Simpson temos∣∣∣∣− h4

180(b− a)f (iv)(η)

∣∣∣∣ ≤ 0.005, i.e.h4

180(1− 0)5e ≤ 0.005



Exercício

Majorante

O majorante de f (iv) no intervalo [0, 1] é então

M = majη∈[0,1]f(iv)(η) = (1 + 4)e1 = 5e

Erro de truncaturaDa fórmula composta de Simpson temos∣∣∣∣− h4

180(b− a)f (iv)(η)

∣∣∣∣ ≤ 0.005, i.e.h4

180(1− 0)5e ≤ 0.005



Exercício

Resolvendo em ordem a h

h4 ≤ 0.005× 1805e

, ⇒ h ≤ 0.507276

Determinar n

Sabendo que h = b−an = 1

n temos

1n≤ 0.507276 ⇔ n ≥ 1.971313

Calculo de n

Como a regra de Simpson só pode ser aplicada para um valor de n partemos então que n = 2 (menor n par que satisfaz a condição) econsequentemente h = 0.5.



Exercício


h4 ≤ 0.005× 1805e

, ⇒ h ≤ 0.507276

Determinar n


n temos

1n≤ 0.507276 ⇔ n ≥ 1.971313

Calculo de n




Exercício


h4 ≤ 0.005× 1805e

, ⇒ h ≤ 0.507276

Determinar n


n temos

1n≤ 0.507276 ⇔ n ≥ 1.971313

Calculo de n




Cálculo de I

Com h = 0.5 temos que

x 0.0 0.5 1.0f(x) 0.000000 0.824361 2.718282

Aplicando a regra de Simpson obtemos:

I ≈ 0.53

(0 + 4× 0.824361 + 2.718282)

= 1.002621


Documents

Métodos Numéricos - norg.uminho.pt · Métodos Numéricos A. Ismael F. Vaz Departamento de Produção e Sistemas Escola de Engenharia Universidade do Minho [email protected] Mestrado