mat lab - explicaçao

Prévia do material em texto

1 
 
Capítulo II – Interpolação Polinomial 
 
 
Manual das Funções: 
 
 Difdiv 
 Horner 
 Interpnewton 
 Polyvalh 
 Vandermond 
 
Função Difdiv 
dd = difdiv(x,y,k) 
 
Dada uma tabela de nós x e valores nodais y, difdiv constrói as diferenças divididas de 
ordem k. 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» dd = difdiv(x,y,k) 
 
x – vetor dos nós, 
y – vetor dos valores nodais, 
k – ordem da diferença dividida. 
 
dd – vetor com as diferenças divididas de ordem k. 
 
Exemplo: 
Dada a tabela: 
x 0.0 2.0 3.0 5.0 6.0 
y 6.0 -1.0 -8.0 -3.0 0.0 
 
Construir a tabela das diferenças divididas até à ordem 3. 
 
2 
 
Sendo: 
x = [0 2 3 5 6] 
y = [6 -1 -8 -3 0] 
k = 1, 2, 3. 
» x = [0 2 3 5 6]; 
» y = [6 -1 -8 -3 0]; 
» dd1 = difdiv(x,y,1) 
 
dd1 = 
 -3.5000 
 -7.0000 
 2.5000 
 3.0000 
 
» dd2 = difdiv(x,y,2) 
 
dd2 = 
 -1.1667 
 3.1667 
 0.1667 
 
» dd3 = difdiv(x,y,3) 
 
dd3 = 
 0.8667 
 -0.7500 
» 
 O que resultaria na seguinte tabela: �� �� [.,.] [.,.,.] [.,.,.,.] 
 
0 6.0 
 -3.5 
2 -1.0 -1.1667 
 7.0 0.8667 
3 -8.0 3.1667 
 2.5 -0.75 
5 -3.0 0.1667 
 3.0 
6 0.0 
 
 
Função Horner 
 
[ncoef, ncentros] = horner(coef,centros,ncentros) 
 
3 
 
 
Dado um polinómio de grau n, na forma de Newton com os seus coeficientes: coef = [��		���	 	…	��] e centros: centros = [
�		
���…	
�], Horner calcula os novos 
coeficientes ncoef = [��� 		���	� 	… 	��� ] do mesmo polinómio com novos centros: 
ncentros = [[
�� 		
���� …	
��], utilizando o algoritmo de Horner. 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
 
» [ncoef, ncentros] = horner(coef,centros,ncent) 
 
 
 coef – vetor dos coeficientes do polinómio na forma de Newton com centros. 
 centros – vetor dos centros do polinómio na forma de Newton com centros. 
 ncentros – vetor do novos centros. 
 
 ncoef – vetor dos coeficientes do polinómio interpolador na forma de Newton com os 
novos centros em ncentros. (coef = [��� 		���	� …	��� ]) 
 ncentros – vetor do novos centros. 
 
 
exemplo: 
 
Dado o polinómio de grau 3 na forma de potências simples: 
���� = −0.25�� + 2.0�� − 3.75� + 1.0 
Com: 
coef = [-0.25 2.0 -3.75 1.0] 
centros = [0 0 0] 
pretende-se apresentar este polinómio com centros em: 
ncentros = [2 -1 5] 
A determinação dos novos coeficientes é obtida da seguinte forma: 
» coef = [-0.25 2.0 -3.75 1]; 
» centros = [0 0 0]; 
» ncentros = [2 -1 5]; 
» [ncoef,ncentros] = horner(coef,centros,ncentros) 
ncoef = 
4 
 
 -0.2500 0.5000 -1.0000 1.0000 
ncentros = 2 -1 5 
» 
 
Onde ncoef é o vetor dos novos coeficientes do polinómio na forma de Newton, o que no 
caso, é representado por: 
 ���� = −0.25�� − 2��� + 1��� − 5� + 0.5�� + 1��� − 5� − 1.0�� − 5� + 1.0 
 
 
Função Interpnewton 
 
 [coef,centros] = interpnewton(x,y) 
 
 
Dada uma tabela com os nós no vetor x e os valores nodais no vetor y, interpnewton 
calcula os centros e os coeficientes do polinómio interpolador na forma de newton com 
centros. 
 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» [coef,centros] = interpnewton(x,y) 
 
 
 x – vetor dos nós; 
 y – vetor dos valores nodais. 
 
coef – vetor dos coeficientes do polinómio interpolador na forma de Newton com centros 
(coef = [��		���	 	…	��]) 
centros – centros do polinómio interpolador na forma de Newton com centros (centros = 
[
�		
���…	
�]) 
 
exemplo: 
 
Dada a tabela: 
 
x 0 1 3 4 
y 1 -1 1 2 
 
A determinação do polinómio interpolador de Newton é obtida da seguinte forma: 
» x = [0 1 3 4]; 
» y = [1 -1 1 2]; 
5 
 
» [coef,centros] = interpnewton(x,y) 
coef = 
 -0.2500 1.0000 -2.0000 1.0000 
centros = 
 3 1 0 
 
Onde coef é o vetor dos coeficientes e centros é o vetor de centros para o polinómio na 
forma de Newton com centros. 
O que no caso, representa o seguinte polinómio: 
 ���� = −0.25�� − 3��� − 1��� − 0� + 1.0�� − 1��� − 0� − 2.0�� − 0� + 1.0 
 
 
Função Polyvalh 
 
y = polyvalh(coef,centros,x) 
 
 
Polyvalh calcula o valor de um polinómio na forma de Newton com centros num ponto x 
dado. 
 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» y = polyvalh(coef,centros,x) 
 
 
coef – vetor dos coeficientes do polinómio; 
centros – vetor dos centros; 
x – valor onde se pretende avaliar o polinómio 
 
y – valor do polinómio na forma de Newton com centros num ponto dado. 
 
exemplo: 
 
Dado o polinómio: 
 ���� = −0.25�� − 3��� − 1��� − 0� + 1.0�� − 1��� − 0� − 2.0�� − 0� + 1.0 
 
Que se pretende avaliar no ponto x = 3.2, ou seja, calcular: 
 ��3.2� 
6 
 
 
Com: 
Coef = [-0.25 1.0 -2.0 1.0] 
centros = [3 1 0] 
 
x = 3.2 
 
o valor de ��3.2� é obtido da seguinte forma: 
» coef = [-0.25 1.0 -2.0 1.0]; 
» centros = [3 1 0]; 
» x = 3.2; 
» y = polyvalh(coef,centros,x) 
y = 
 1.28800 
» 
Sendo ��3.2� =1.288. 
 
Função Vandermond 
 
coef = vandermond(x,y) 
 
 
Dada uma tabela com os nós no vetor x e os valores nodais no vetor y, vandermond 
calcula os coeficientes do polinómio interpolador na forma de potências simples 
recorrendo ao sistema de Vandermond. 
 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
 
» coef = vandermond(x,y) 
 
 
 x – vetor dos nós; 
 y – vetor dos valores nodais. 
 
7 
 
coef – vetor dos coeficientes do polinómio interpolador na forma de potências simples 
(coef = [��		���	 	…	��]) 
 
exemplo: 
 
Dada a tabela: 
 
x 0 1 3 4 
y 1 -1 1 2 
 
A determinação do polinómio interpolador de Newton é obtida da seguinte forma: 
» x = [0 1 3 4]; 
» y = [1 -1 1 2]; 
» coef = vandermond(x,y) 
coef = 
 -0.2500 2.0000 -3.7500 1.0000 
 
Onde coef é o vetor dos coeficientes do polinómio na forma de potências simples. 
O que no caso, representa o seguinte polinómio: 
 ���� = −0.25�� + 2.0�� − 3.75� + 1.0 
 
 
 
 
Capítulo III– Solução de equações não lineares 
 
 
Manual das Funções: 
 
 Bisseccao 
 Newton 
 Secante 
 
 
Função Bisseccao 
xk = bisseccao(a,b,tol) 
8 
 
Dados os limites a e b onde existe uma raiz da função f(x), a definir, bisseccao determina 
uma aproximação da raiz, pelo método da bissecção, nesse intervalo com o majorante do 
erro absoluto tol especificado. 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» y = bisseccao(a,b,tol) 
 
a – limite inferior do intervalo que contem a raiz, 
b – limite superior do intervalo que contem a raiz, 
tol – majorante do erro. 
 
xk – valor aproximado da raiz. 
 
Exemplo: 
Dada a função: 
!��� = �" + 2� + 1.0 
que se sabe ter um zero no intervalo [-1.5, -1], determinar uma aproximação dessa raiz 
pelo método da bissecção, com um majorante do erro de 0.5 × 10��. 
Primeiro é necessário definir a função !��� dada. Isso é feito acedendo ao ficheiro 
bisseccao.m e no corpo da função f(x) já existente no final do ficheiro, introduzindo a 
linha: 
 y=x^4+2*x+1; 
 
tal que o resultado seja o seguinte: 
 
 
function y = f(x) 
 
 y = x^4 + 2*x + 1; 
 
end 
Alterado assim o ficheiro bisseccao.m, deve proceder-se à sua gravação em disco. 
 
Sendo: 
a = -1.5 
b = -1.0 
9 
 
tol =	0.5 × 10�� 
O valor aproximado é obtido com a introdução dos seguintes dados, e apóscada iteração 
deve premir-se uma tecla qualquer para que sejam apresentados os dados da próxima 
iteração: 
» a = -1.5 
» b = -1.0; 
» tol = 0.5e-2; 
» y = bisseccao(a,b,tol) 
 k ak bk xk f(ak) f(xk) f(ak)*f(xk) 
 0 -1.50000000 -1.00000000 -1.25000000 1.06 -1.06e+000 <0 
 1 -1.50000000 -1.25000000 -1.37500000 1.06 -1.76e-001 <0 
 2 -1.50000000 -1.37500000 -1.43750000 1.06 3.95e-001 >0 
 3 -1.43750000 -1.37500000 -1.40625000 0.40 9.82e-002 >0 
 4 -1.40625000 -1.37500000 -1.39062500 0.10 -4.15e-002 <0 
 5 -1.40625000 -1.39062500 -1.39843750 0.10 2.76e-002 >0 
 6 -1.39843750 -1.39062500 -1.39453125 
 
y = 
 -1.39453125 
Sendo o valor aproximado da raiz, dado por: 
 �$ =-1.394 
Que se sabe ter 2 casas decimais significativas. 
 
 
Função Newton 
y = Newton(x0,tol,maxiter) 
Dada uma estimativa inicial x0 da raiz da função f(x) a definir, e os critérios de paragem 
tol e maxiter, newton, determina uma aproximação dessa raiz. tol é o valor, tal que 
quando o valor da diferença, em módulo, entre duas iterações consecutivas seja inferior 
ou igual, o método termina. maxiter é o número máximo de iterações a partir do qual o 
método é terminado. 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» y = newton(x0,tol,maxiter) 
 
10 
 
x0 – estimativa inicial da raiz, 
tol – critério de paragem, tal que quando se verifique que a diferença, em módulo, entre 
duas iterações consecutivas não excede esse valor: |�&'� − �&| ≤ )*+,o método termina. 
maxiter – número máximo de iterações que o método deve realizar, caso não aconteça 
que durante esse número de iterações, se tenha verificado: |�&'� − �&| ≤ )*+. 
 
y – valor aproximado da raiz. 
 
Exemplo: 
Dada a função: 
!��� = ,- + �� − 2.0 
que se sabe ter um zero no intervalo [-1.4, -1.3], determinar uma aproximação dessa raiz 
pelo método de Newton, com um critério de paragem tal que |�&'� − �&| ≤ 0.5 × 10�$. 
Primeiro é necessário definir a função !��� dada e a sua derivada, sendo !′��� = ,- +2�. 
Isso é feito acedendo ao ficheiro newton.m e no corpo das função f(x) e df(x) já 
existentes no final do ficheiro, introduzindo as linhas: 
 yf=exp(x)+x^2*x-2; 
 e 
 ydf=exp(x)+2*x; 
 
tal que o resultado seja o seguinte: 
 
 
function yf = f(x) 
 
 yf = exp(x) + x^2 - 2; 
 
end 
 
 
function ydf = df(x) 
 
 ydf = exp(x) + 2*x; 
 
end 
 
Alterado assim o ficheiro newton.m, deve proceder-se à sua gravação em disco. 
 
Sendo: 
11 
 
x0 = -1.4 (ou qualquer outro valor no intervalo [-1.4, -1.3]). 
tol = 0.5 × 10�$ 
e atribuindo, por exemplo, um valor de 100 para o número máximo de 
iterações 
maxiter=100 
 
O valor aproximado é obtido com a introdução dos seguintes dados: 
» x0 = -1.4 
» tol = 0.5e-6; 
» maxiter = 100; 
» y = newton(x0,tol,maxiter) 
 
x1= -1.3999999999999999 + 0.0809104403120486 = -1.3190895596879513 
x2= -1.3190895596879513 + 0.0031111390233119 = -1.3159784206646394 
x3= -1.3159784206646394 + 0.0000046428580067 = -1.3159737778066327 
x4= -1.3159737778066327 + 0.0000000000103426 = -1.3159737777962901 
 y = 
 -1.315973777796290 
 
O que se pode concluir, neste caso, observando a sucessão de valores obtidos, que o 
valor �" = −1.315973777 (com todos os dígitos significativos) é uma aproximação 
da raiz com pelo menos 9 casas decimais significativas. Isto porque, se observarmos, 
de �� para �", as primeiras 9 casas decimais de �� já não se alteraram, logo �� tem 
essas 9 casas decimais significativas, o que implica que �" tem pelo menos aquelas 9 
casas decimais significativas. Ou seja, pode-se concluir à priori que o erro associado a �" é inferior a 0.5 × 10�0, isto porque a sucessão de valores gerada pelo método está 
a ser convergente para a raiz. 
 
 
Função Secante 
y = secante(x0,x,tol,maxiter) 
 
12 
 
Dada duas estimativas iniciais x0 e x, com x0 ≠ x, a raiz da função f(x) a definir, e os 
critérios de paragem tol e maxiter, secante, determina uma aproximação dessa raiz. tol é 
o valor, tal que quando o valor da diferença, em módulo, entre duas iterações 
consecutivas seja inferior ou igual, o método termina. maxiter é o número máximo de 
iterações a partir do qual o método é terminado. 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» y = secante(x0,x,tol,maxiter) 
 
x0 – estimativa inicial da raiz, 
x – estimativa inicial da raiz, 
tol – critério de paragem, tal que quando se verifique que a diferença, em módulo, entre 
duas iterações consecutivas não excede esse valor: |�&'� − �&| ≤ )*+,o método termina. 
maxiter – número máximo de iterações que o método deve realizar, caso não aconteça 
que durante esse número de iterações, se tenha verificado: |�&'� − �&| ≤ )*+. 
 
y – valor aproximado da raiz. 
 
Exemplo: 
Dada a função: 
!��� = ,- + �� − 2.0 
que se sabe ter um zero no intervalo [-1.4, -1.3], determinar uma aproximação dessa raiz 
pelo método da Secante, com um critério de paragem tal que |�&'� − �&| ≤ 0.5 × 10�$. 
Primeiro é necessário definir a função !��� dada. 
Isso é feito acedendo ao ficheiro secante.m e no corpo da função f(x) já existente no final 
do ficheiro, introduzindo a linha: 
 y=exp(x)+x^2*x-2; 
 
 
tal que o resultado seja o seguinte: 
 
 
function y = f(x) 
 
 y = exp(x) + x^2 - 2; 
 
end 
Alterado assim o ficheiro secante.m, deve proceder-se à sua gravação em disco. 
13 
 
 
Sendo: 
x0 = -1.4 (ou qualquer outro valor no intervalo [-1.4, -1.3]). 
x = -1.3 (ou qualquer outro valor no intervalo [-1.4, -1.3] diferente de 
x0) 
tol = 0.5 × 10�$ 
e atribuindo, por exemplo, um valor de 100 para o número máximo de 
iterações 
maxiter=100 
O valor aproximado é obtido com a introdução dos seguintes dados: 
» x0 = -1.4; 
» x0 = -1.3; 
» tol = 0.5e-6; 
» maxiter = 100; 
» y = secante(x0,x,tol,maxiter) 
 
x1= -1.3000000000000000 - 0.0153517221759483 = -1.3153517221759483 
x2= -1.3153517221759483 - 0.0006268645817516 = -1.3159785867576999 
x3= -1.3159785867576999 + 0.0000048103971465 = -1.3159737763605535 
x4= -1.3159737763605535 - 0.0000000014357334 = -1.3159737777962868 
y = 
 -1.315973777796287 
 
O que se pode concluir, neste caso, observando a sucessão de valores obtidos, que o 
valor �" = −1.31597377 (com todos os dígitos significativos) é uma aproximação 
da raiz com pelo menos 8 casas decimais significativas. Isto porque, se observarmos, 
de �� para �", as primeiras 8 casas decimais de �� já não se alteraram, logo �� tem 
essas 8 casas decimais significativas, o que implica que �" tem pelo menos aquelas 8 
casas decimais significativas. Ou seja, pode-se concluir à priori que o erro associado a �" é inferior a 0.5 × 10�2, isto porque a sucessão de valores gerada pelo método está 
a ser convergente para a raiz. 
 
14 
 
Capítulo IV– Integração numérica 
 
 
Manual das Funções: 
 
 newton_cotes 
 i_gauss 
 
 
Função Newton_cotes 
y = newton_cotes(a,b,n,m,fx) 
Dada uma função f(x) na forma analítica ou na forma de vetor fx, newton_cotes calcula o 
valor aproximado do integral definido, da função no intervalo [a,b],utilizando as regras 
de Newton-Cotes com um polinómio interpolador de grau n, com m intervalos. 
 
A chamada à função, em ambiente Matlab ou Octave, pode ser feita de duas formas, 
digitando na linha de comandos quando a função é definida na forma analítica: 
 
» y = newton_cotes(a,b,n,m) 
 
 
ou, quando a função é definida na forma de vetor:» y = newton_cotes(a,b,n,m,fx) 
 
 
a – limite inferior de integração, 
b – limite superior de integração, 
n – grau do polinómio interpolador, ou número de nós + 1, equidistantes no intervalo de 
integração [a, b], 
m – número de intervalos, 
fx – vetor contendo os valores da função a integrar. 
 
y – valor aproximado do integral definido. 
 
Exemplo: 
Dada a função: 
!��� = ��,345	�-� 
15 
 
Pretende-se determinar um valor aproximado de 6 !���7��8 utilizando as regras de 
Newton cotes com n=2: 
a) Simples, 
b) Composta, com 5 intervalos. 
Primeiro é necessário definir a função !��� dada. 
Isso é feito acedendo ao ficheiro newton_cotes.m e no corpo da função f(x) já existente 
no final do ficheiro, introduz-se a linha: 
 y= x^2*exp(cos(x)); 
 
 
tal que o resultado seja o seguinte: 
 
 
function y = f(x) 
 
 y = x^2 * exp(cos(x)); 
 
end 
Alterado assim o ficheiro newton_cotes.m, deve proceder-se à sua gravação em disco. 
a) 
Sendo: 
a = 0, b = 1, n = 2 e m = 1 
O valor aproximado é obtido com a introdução dos seguintes dados: 
» a = 0; b = 2.0; n = 2; m = 1; 
» y = newton_cotes(a,b,n,m) 
y = 
 3.1681 
Em que y=3.1681, corresponde ao valor aproximado do integral definido, obtido pelas 
regras de Simpson simples. 
 
b) neste caso, apenas o valor de m é alterado para m = 5; 
O valor aproximado é obtido com a introdução dos seguintes dados: 
» a = 0; b = 2.0; n = 2; m = 5; 
16 
 
» y = newton_cotes(a,b,n,m) 
y = 
 3.021308473305821 
Em que y=3.021, representa ao valor aproximado do integral definido, obtido pelas regras 
de Simpson composta com 5 intervalos. 
Para sabermos o número de casas decimais significativas, poderíamos aplicar as regras de 
Newton Cotes, com um número de intervalos suficientemente grande e tirar daí 
conclusões: 
» a = 0; b = 2.0; n = 2; 
» y=newton_cotes(a,b,n,5) 
y = 
 3.021308473305821 
» y=newton_cotes(a,b,n,10) 
y = 
 3.021277881251750 
» y=newton_cotes(a,b,n,20) 
y = 
 3.021275973373415 
» y=newton_cotes(a,b,n,30) 
y = 
 3.021275871369385 
Pode-se afirmar, pela observação da sucessão de valores que o resultado pedido, tem na 
alínea a) 0 c.d.s. e na alínea b) 3 c.d.s. 
 
Caso se pretenda calcular 6 !���7�9: em que f(x) é dada na forma de tabela: 
x 0 1 2 3 4 5 6 
f(x) 1 -1 1 2 3 7 5 
 
Considerando que se pretende a regra de Newton Cotes composta com n=2 e 
m=3. 
 
Sendo a=0, b=6 e fx=[1 -1 1 2 3 7 5], vem: 
 
 
» a = 0; b = 6.0; n = 2; m = 3; 
» fx=[1 -1 1 2 3 7 5]; 
» y=newton_cotes(a,b,n,m,fx) 
y = 
 
 15.333333333333334 
 
Com o valor aproximado do integral dado por y=15.333. 
 
17 
 
Quanto ao número de c.d.s. para o resultado, não é de todo possível 
determinar, sem que seja dado uma informação acerca da função. 
De notar que o valor aproximado do integral assenta na disposição 
equidistante dos nós. 
 
 
Função i_gauss 
y = i_gauss(n) 
Dada uma função f(x) na forma analítica, i_gauss calcula o valor aproximado do integral 
definido, da função no intervalo normalizado [-1,1],utilizando a quadratura de Gauss 
Legendre com n nós. 
 
 
A chamada à função, em ambiente Matlab ou Octave, pode ser feita de duas formas, 
digitando na linha de comandos: 
 
» y = i_gauss(n) 
 
 
n – número de nós de integração, que são os n zeros dos polinómios de Legendre de grau 
n. 
 
y – valor aproximado do integral definido. 
 
Exemplo: 
Dada a função: 
!��� = ��,345	�-� 
Pretende-se determinar um valor aproximado de 6 !���7���� utilizando as regras de 
Gauss Legendre com n=4: 
Primeiro é necessário definir a função !��� dada. 
Isso é feito acedendo ao ficheiro i_gauss.m e no corpo da função f(x) já existente no final 
do ficheiro, introduzir-se a linha: 
 y= x^2*exp(cos(x)); 
 
 
tal que o resultado seja o seguinte: 
18 
 
 
 
function y = f(x) 
 
 y = x^2 * exp(cos(x)); 
 
end 
Alterado assim o ficheiro i_gauss.m, deve proceder-se à sua gravação em disco. 
n = 4 
O valor aproximado é obtido com a introdução dos seguintes dados: 
» n = 4; 
» y = i_gauss(n) 
y = 
 1.376808182585136 
Em que y=1.3768 é o resultado obtido. 
Para sabermos o número de casas decimais significativas, poderíamos aplicar a regra 
de Gauss com n seguinte e tirar daí conclusões: 
» y= i_gauss(4) 
 
y = 
 
 1.376808182585136 
 
» y= i_gauss(5) 
y = 
 1.375909787165000 
» y= i_gauss(6) 
y = 
 1.375957137393892 
Pode-se afirmar, pela observação da sucessão de valores que o resultado pedido y= 
1.37 tem 2 c.d.s. 
 
 
19 
 
Capítulo V – Método dos Mínimos Quadrados 
 
Manual das Funções: 
 
 Min_quad 
 
 
 
Função min_quad 
[G,a] = min_quad(x,y) 
Para uma função !��� dada na forma discreta sobre um conjunto de valores expresso no 
vetor x, a função min_quad determina os parâmetros da função aproximadora ;��� dada 
por: 
 
;��� = <=&>&���?&@8 
 
que melhor se ajusta à função !���, no sentido dos mínimos quadrados. 
 
De uma forma geral, o conjunto dos parâmetros =&, k=1…A, dado no vetor a, é 
encontrado resolvendo o sistema de equações normais: 
 B × � = C 
 
em que: 
 B = DE × D 
 
e 
 C = DE × � 
 
sendo: 
 D = [>8��� >���� … >?���] 
 
onde >&��� são as funções base. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» [G,a] = min_quad(x,y) 
 
 
20 
 
x – vetor dos nós, onde a função toma valores, 
y – vetor dos valores nodais da função a aproximar, 
 
G – matriz do sistema normal 
a – vetor dos parâmetros 
 
Exemplo: 
Dada a função a aproximar que se apresenta na tabela: 
 
x 0 1 2 3 4 
y 1 -1 0.5 1 2 
 
e cujo gráfico se apresenta: 
 
Pretende-se determinar a reta como função aproximadora que melhor se ajusta ao 
conjunto de pontos. Ou seja, determinar os parâmetros =8 e =� da função: 
;��� = <=&>&��� =�&@8 =8 + =�	� 
Com as funções base >8��� = 1 e >���� = � 
É necessário aceder ao ficheiro min_quad.m para definir as duas funções de base >8��� e >���� e proceder à construção da matriz X. 
Isso é conseguido, introduzindo as linhas: 
21 
 
 fi0 = x.^0; 
 fi1 = x.^1; 
 
 X = [fi0 fi1]; 
 
tal que o resultado seja o seguinte: 
 
 
% entrar aqui as funcoes base 
 
 fi0 = x.^0; 
 fi1 = x.^1; 
 
 X = [fi0 fi1]; 
 
% terminar aqui 
 
 
 Alterado assim o ficheiro min_quad.m, deve proceder-se à sua gravação em disco. 
De notar que as funções base >&��� têm de ser avaliadas elemento a elemento, daí que 
as operações tenham de ser precedidas de um “.”, refletindo isso mesmo. 
A matriz do sistema e os parâmetros são obtidos com a introdução dos seguintes 
dados: 
» x=[0 1 2 3 4]; y=[1 -1 0.5 1 2]; 
» [G,a] = min_quad(x,y) 
 
G = 
 
 5 10 
 10 30 
 
 
a = 
 
 -0.1000 0.4000 
O vetor de parâmetros a é constituído neste caso pelos parâmetros =8 e =� com os valores 
-0.1 e 0.4 respetivamente, constituindo a reta aproximadora: 
;��� = −0.1 + 0.4	� 
 
 O gráfico mostrando o ajustamento obtido pode ser agora visualizado: 
22 
 
 
Seguem-se exemplos de outras funções aproximadoras e sua implementação: 
Ex. 2 
;��� = <=&>&��� =�&@8 =8 + =�� + =�	�� 
 
 fi0 = x.^0; 
 fi1 = x.^1; 
 fi2 = x.^2; 
 
 X = [fi0 fi1 fi2]; 
 
tal que o resultado seja o seguinte: 
 
 
% entrar aqui as funcoes base 
 
 fi0 = x.^0; 
 fi1 = x.^1; 
 fi2 = x.^2; 
 
 X = [fi0 fi1 fi2]; 
 
% terminar aqui 
Ex. 3 
;��� = <=&>&��� =�&@8 =81 + � + =�,- + =�	sin�	��� 
23 
 
 
 fi0 = 1./(1 + x);fi1 = exp(x); 
 fi2 = sin(x).^2; 
 
 X = [fi0 fi1 fi2]; 
 
tal que o resultado seja o seguinte: 
 
 
% entrar aqui as funcoes base 
 
 fi0 = 1./(1+x); 
 fi1 = exp(x); 
 fi2 = sin(x).^2; 
 
 X = [fi0 fi1 fi2]; 
 
% terminar aqui 
 
24 
 
Capítulo VI – Solução de Sistemas de Equações 
 
 
 
JKLMNOPL	QN	RSTPçõNL	
WXX
YX
XZ[K\NP]NL		 ^ Métodos	diretos�f=)*ghi=çã*	kl�	m
*n7,no=çã*	7,	p=qooré)*7*o	s*A�=s)*o	 tu**+h)+,
g*q)Métodos	iterativos	 tx=s*yhp=qoo − z,h7,+
{ã|	[K\NP]NL − Método	de	Newton
 
 
 
A solução de sistemas de equações lineares através dos métodos diretos, recorre à 
factorização LU da Matriz A do sistema de equações Ax=b, tal que A=L×U. Após a 
factorização da matriz A, a solução do sistema x é obtida pela resolução de dois sistemas 
triangulares, tal que: 
 



=
=
yxU
byL
 
 
A factorização LU pode ser obtida pela condensação de Gauss ou pelos métodos 
compactos de Doolitle ou de Crout. 
Uma vez obtida a factorização LU da matriz A, a solução do sistema Ax=b, para 
diferentes vectores b, pode ser obtida pela resolução dos dois sistemas triangulares, para 
cada vetor b, onde o processo da factorização LU é realizado apenas uma vez, o que 
acarreta redução nos cálculos, pois este processo da factorização LU é o mais pesado em 
termos computacionais, com o número de operações aritméticas da ordem de n� face a n� 
para a resolução dos sistemas triangulares, onde n a ordem do sistema. 
 
Com o objectivo de reduzir os erros de arredondamento, convém escolher para candidatos 
a elementos pivô, elementos com o maior valor absoluto. Isto é conseguido com duas 
estratégias de procura de candidatos a elemento pivô: 
 
 - Pivô parcial 
 - Pivô total 
 
Na estratégia do pivô parcial, apenas se recorre à troca de linhas, e na estratégia de pivô 
total recorre-se à troca de linhas e ou colunas, onde a troca de linhas é expressa na matriz 
de permutação de linhas P e a troca de colunas expressa na matriz de permutação de 
colunas Q. 
 
Assim, quando se utilizou a estratégia de pivô parcial a factorização LU para a obtenção 
das matrizes L e U, terá de verificar-se que L×U=P×A. Para encontrar a solução x do 
sistema, terá que se realizar a resolução de dois sistemas triangulares de acordo com: 
 
25 
 
L y Pb
U x y
=
=



 
 
Quando se utilizou a estratégia de pivô total, a factorização LU para a obtenção das 
matrizes L e U, terá de verificar-se agora que L×U=P×A×Q. Para encontrar a solução x 
do sistema, terá que se realizar a resolução de dois sistemas triangulares e a multiplicação 
de uma matriz por um vetor, de acordo com: 
 
L z Pb
U y z
x Q y
=
=
=





 
 
Uma forma atraente de resolver um sistema de equações lineares, pode ser obtido pelos 
métodos iterativos de Jacobi e Gauss-Seidel. Isto tem grande interesse quando a matriz 
do sistema é de diagonal estritamente dominante por linhas. 
 
Para a obtenção da solução de sistemas de equações não lineares, um método atraente é o 
método de Newton utilizado na solução de equações não lineares, agora adaptado, onde 
se recorre à matriz Jacobiana J(x). O método iterativo de Newton para sistemas de 
equações não lineares pode ser realizado de acordo com: 
 
( ) ( )



∆+=
−=∆
+ )()()1(
)()()(
kkk
kkk
xxx
xfxxJ
 
 
Especificando como critério de paragem: 
 
ε≤−+
p
kk )()1( xx 
 
 
1 – Sistemas de Equações Lineares 
 
a) – Métodos diretos 
 
axb 
axb_pp 
axb_pt 
lu_crout 
lu_doolitle 
lu_gauss 
lu_pp 
lu_pt 
lxb 
uxb 
 
26 
 
 
b) – Métodos iterativos 
g_seidel 
jacobi 
 
2 – Sistemas de Equações Não Lineares 
 
 soleqnl_inter 
 soleqnl 
 
 
1 – Sistemas de Equações não lineares – Métodos diretos 
 
Considere-se o seguinte sistema de equações lineares Ax=b: 
 
 1 2 −12 3 0−2 1 3 € 
������€ = 
−3−28 € 
 
 
Função Axb 
x = axb(A,b) 
 
Dada um sistema de equações lineares Ax=b, axb calcula o vetor solução x, dada a matriz 
do sistema A e o vetor dos termos independentes b, pela fatorização LU obtida pela 
condensação de Gauss e solução dos sistemas triangulares inferior e superior. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» x = axb(A,b) 
 
 
 
A – Matriz do sistema, 
b – Vetor dos termos independentes, 
 
x – vetor solução. 
 
Exemplo: 
Para o sistema dado, cuja solução se pretende determinar: 
A solução é obtida com a introdução dos seguintes dados: 
27 
 
» A=[1 2 -1; 2 3 0; -2 1 3]; b=[-3 -2 8]’; 
» x = axb(A,b) 
x = 
 -1 0 2 
Em que ������€ = 
−102 € corresponde à solução do sistema. 
 
Função Axb_pp 
x = axb_pp(A,b) 
 
Dada um sistema de equações lineares Ax=b, axb calcula o vetor solução x, dada a matriz 
do sistema A e o vetor dos termos independentes b, pela fatorização LU obtida pela 
condensação de Gauss e estratégia de pivô parcial. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» x = axb_pp(A,b) 
 
 
 
A – Matriz do sistema, 
b – Vetor dos termos independentes, 
 
x – vetor solução. 
 
Exemplo: 
Para o sistema dado, cuja solução se pretende determinar: 
A solução é obtida com a introdução dos seguintes dados: 
» A=[1 2 -1; 2 3 0; -2 1 3]; b=[-3 -2 8]’; 
» x = axb_pp(A,b) 
x = 
 -1 0 2 
28 
 
Em que ������€ = 
−102 € corresponde à solução do sistema. 
 
Função Axb_pt 
x = axb_pt(A,b) 
 
Dada um sistema de equações lineares Ax=b, axb calcula o vetor solução x, dada a matriz 
do sistema A e o vetor dos termos independentes b, pela fatorização LU obtida pela 
condensação de Gauss e estratégia de pivô total. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» x = axb_pt(A,b) 
 
 
 
A – Matriz do sistema, 
b – Vetor dos termos independentes, 
 
x – vetor solução. 
 
Exemplo: 
Para o sistema dado, cuja solução se pretende determinar: 
A solução é obtida com a introdução dos seguintes dados: 
» A=[1 2 -1; 2 3 0; -2 1 3]; b=[-3 -2 8]’; 
» x = axb_pt(A,b) 
x = 
 -1 0 2 
Em que ������€ = 
−102 € corresponde à solução do sistema. 
 
 
29 
 
Função Lu_crout 
[L,U] = lu_crout(A) 
 
Dada a matriz A do sistema lu_crout realiza a factorização LU pelo método compacto de 
Crout . 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» [L,U] = lu_crout(A) 
 
 
 
A – Matriz do sistema, 
 
 
L – Matriz triangular inferior, 
U – Matriz triangular superior, 
 
Exemplo: 
Considere-se o sistema dado, cuja factorização LU da matriz A do sistema se pretende 
obter pelo método compacto de Crout: 
A resultado é obtido com a introdução dos seguintes dados: 
» A=[1 2 -1; 2 3 0; -2 1 3];b=[-3 -2 8]’; 
» [L,U]=lu_crout(A) 
L = 
 1 0 0 
 2 -1 0 
 -2 5 11 
U = 
 1 2 -1 
 0 1 -2 
 0 0 1 
 
Em que k =  1 0 02 −1 0−2 5 11€ 	e	 l = 
1 2 −10 1 −20 0 1 € é a factorização LU obtida pelo 
método compacto de Crout. Repare-se que a diagonal da matriz U é constituída por 
valores unitários. 
30 
 
 
Função Lu_doolitle 
[L,U] = lu_doolitle(A) 
 
Dada a matriz A do sistema, lu_doolitle realiza a factorização LU pelo método compacto 
de Doolitle. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» [L,U] = lu_doolitle(A) 
 
 
A – Matriz do sistema, 
 
 
L – Matriz triangular inferior, 
U – Matriz triangular superior, 
 
Exemplo: 
Dado o sistemade equações acima, cuja factorização LU da matriz A se pretende obter 
pelo método compacto de Doolitle. 
O resultado é obtido com a introdução dos seguintes dados: 
» A=[1 2 -1; 2 3 0; -2 1 3]; 
 
» [L,U]=lu_doolitle(A) 
 
L = 
 
 1 0 0 
 2 1 0 
 -2 -5 1 
 
 
U = 
 
 1 2 -1 
 0 -1 2 
 0 0 11 
 
 
31 
 
Em que k =  1 0 02 1 0−2 −5 1€ 	e	 l = 
1 2 −10 −1 20 0 11€ é a factorização LU obtida pelo 
método compacto de Doolitle. Repare-se que é agora a diagonal da matriz L que é 
constituída por valores unitários. 
 
Função Lu_gauss 
[L,U] = lu_gauss(A) 
 
Dada a matriz A do sistema, lu_gauss realiza a factorização LU pelo método de Gauss. 
 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» [L,U] = lu_gauss(A) 
 
 
A – Matriz do sistema, 
 
 
L – Matriz triangular inferior, 
U – Matriz triangular superior, 
 
Exemplo: 
Dado o sistema de equações, cuja factorização LU da matriz A se pretende obter pela 
condensação de Gauss. 
A resultado é obtido com a introdução dos seguintes dados: 
 
 
» A=[1 2 -1; 2 3 0; -2 1 3]; 
» [L,U,x]=lu_gauss(A,b) 
 
L = 
 
 1 0 0 
 2 1 0 
 -2 -5 1 
 
32 
 
 
U = 
 
 1 2 -1 
 0 -1 2 
 0 0 11 
 
 
Em que k =  1 0 02 1 0−2 −5 1€ 	e	 l = 
1 2 −10 −1 20 0 11€ é a factorização LU obtida pelo 
método de Gauss. Repare-se que é agora a diagonal da matriz L que é constituída por 
valores unitários, em semelhança ao método de Doolitle. 
 
Função Lu_gauss_inter 
LU = lu_gauss_inter(A) 
 
Dada a matriz A do sistema, lu_gauss_inter realiza a factorização LU pelo método de 
Gauss de uma forma interativa, podendo-se assim visualizar os passos intermédios do 
processo de condensação da matriz A, para cada sub-matriz ativa, podendo-se visualizar 
os multiplicadores que vão construindo a matriz L triangular inferior. A matriz final tem 
concatenado as duas matrizes triangulares L e U, onde os elementos da diagonal de L são 
unitários. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» LU = lu_gauss_inter(A) 
 
 
A – Matriz do sistema, 
 
 
LU – Matriz que concatena as matrizes triangulares Le U 
 
Exemplo: 
Dada a matriz do sistema acima, cuja factorização LU se pretende obter de uma forma 
interativa. 
A resultado é obtido com a introdução dos seguintes dados: 
» A=[1 2 -1; 2 3 0; -2 1 3]; 
33 
 
» LU=lu_gauss_inter(A) 
 
a = 
 
 1 2 -1 
 2 -1 2 
 -2 5 1 
 
 
a = 
 
 1 2 -1 
 2 -1 2 
 -2 -5 11 
 
 
LU = 
 
 1 2 -1 
 2 -1 2 
 -2 -5 11 
 É visível os dois passos intermédios obtidos durante a condensação e a construção da 
matriz triangular inferior L através dos multiplicadores. 
A interpretação final do resultado tem em conta que a matriz LU tem em si concatenadas 
as duas matrizes triangulares L e U, em que k =  1 0 02 1 0−2 −5 1€ 	e	 l = 
1 2 −10 −1 20 0 11€ é 
a factorização LU obtida pelo método de Gauss. Repare-se que é agora a diagonal da 
matriz L que é constituída por valores unitários, em semelhança ao método de Doolitle. 
 
Função Lu_pp 
[L,U,p] = lu_pp(A) 
 
Dada a matriz A do sistema, lu_pp realiza a factorização LU pelo método de Gauss com 
estratégia de pivô parcial. 
 
O vetor p é o vetor permutação de linhas. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» [L,U,p] = lu_pp(A) 
 
 
A – Matriz do sistema, 
 
 
L – Matriz triangular inferior, 
34 
 
U – Matriz triangular superior, 
p – vetor de permutação de linhas, 
 
Exemplo: 
Dado o sistema de equações acima, cuja factorização LU com estratégia de pivô parcial 
da matriz A se pretende obter: 
A resultado é obtido com a introdução dos seguintes dados: 
» A=[1 2 -1; 2 3 0; -2 1 3]; 
» [L,U,p]=lu_pp(A) 
 
L = 
 
 1.0000 0 0 
 -1.0000 1.0000 0 
 0.5000 0.1250 1.0000 
 
 
U = 
 
 2.0000 3.0000 0 
 0 4.0000 3.0000 
 0 0 -1.3750 
 
 
p = 
 
 2 
 3 
 1 
 
 
Em que k =  1 0 0−1 1 00.5 0.125 1€ 	e	 l = 
2 3 −10 4 30 0 −1.375€ é a factorização LU obtida pelo 
método de Gauss com estratégia de pivô parcial, � = 231€ é o vetor de permutação de 
linhas. Repare-se que é a diagonal da matriz L que é constituída por valores unitários. 
De notar que o produto de L por U irá resultar na matriz inicial A com as linhas trocadas 
de acordo com a indexação referida no vetor p de permutação de linhas, tal como se pode 
visualizar: 
» L*U 
 
ans = 
 
 2 3 0 
 -2 1 3 
 1 2 -1 
35 
 
 
» A(p,:) 
 
ans = 
 
 2 3 0 
 -2 1 3 
 1 2 -1 
 
Função Lu_pt 
[L,U,p,q] = lu_pp(A) 
 
Dada a matriz A do sistema, lu_pt realiza a factorização LU pelo método de Gauss com 
estratégia de pivô total. 
 
Em que p é o vetor permutação de linhas e q é o vetor permutação de colunas. 
 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» [L,U,p,q,x] = lu_pt(A,b) 
 
 
A – Matriz do sistema, 
 
 
L – Matriz triangular inferior, 
U – Matriz triangular superior, 
p – vetor de permutação de linhas, 
q – vetor de permutação de colunas, 
 
Exemplo: 
Dado o sistema de equações acima, cuja factorização LU com estratégia de pivô total da 
matriz A se pretende obter: 
A resultado é obtido com a introdução dos seguintes dados: 
» A=[1 2 -1; 2 3 0; -2 1 3]; 
 
» [L,U,p,q]=lu_pt(A) 
 
L = 
 
 1.0000 0 0 
36 
 
 0.3333 1.0000 0 
 0.6667 -0.3333 1.0000 
 
 
U = 
 
 3.0000 0 2.0000 
 0 3.0000 -2.6667 
 0 0 -1.2222 
 
 
p = 
 
 2 
 3 
 1 
 
 
q = 
 
 2 3 1 
 
 
 
Em que k =  1 0 00.3333 1 00.6667 −0.3333 1€ 	e	 l = 
3 0 20 3 −2.66670 0 −1.2222€ é a factorização LU 
obtida pelo método de Gauss com estratégia de pivô total, � = 231€ é o vetor de 
permutação de linhas, ƒ = [2 3 1]. 
Repare-se que é a diagonal da matriz L que é constituída por valores unitários. 
De notar que o produto de L por U irá resultar na matriz inicial A, mas agora com as 
linhas e colunas trocadas de acordo com as indexações referidas nos vetores p e q de 
permutação de linhas e colunas, respectivamente, tal como se pode visualizar: 
» L*U 
 
ans = 
 
 3.0000 0 2.0000 
 1.0000 3.0000 -2.0000 
 2.0000 -1.0000 1.0000 
 
» A(p,q) 
 
ans = 
 
 3 0 2 
 1 3 -2 
 2 -1 1 
 
37 
 
Função lxb 
x = lxb(L,b) 
 
Dada um sistema de equações triangular inferior Lx=b, lxb calcula o vetor solução x, 
pelo método das substituições descendentes. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» x = lxb(L,b) 
 
 
 
L – Matriz do sistema triangular inferior, 
b – Vetor dos termos independentes, 
 
x – vetor solução. 
 
Exemplo: 
Dado o sistema triangular inferior: 
 
 1 0 02 1 0−2 −5 1€ 
������€ = 
−3−28 € 
 
A solução é obtida com a introdução dos seguintes dados: 
» L=[1 0 0; 2 1 0; -2 -5 3]; b=[-3 -2 8]’; 
» x = lxb(L,b) 
x = 
 -3 4 22 
Em que ������€ = 
−3422€ corresponde à solução do sistema. 
 
Função uxb 
38 
 
x = uxb(U,b)Dada um sistema de equações triangular superior Ux=b, uxb calcula o vetor solução x, 
pelo método das substituições ascendentes. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» x = uxb(U,b) 
 
 
 
U – Matriz do sistema triangular inferior, 
b – Vetor dos termos independentes, 
 
x – vetor solução. 
 
Exemplo: 
Dado o sistema triangular superior: 
 
1 2 −10 −1 20 0 11€ 
������€ = 
−3422€ 
 
A solução é obtida com a introdução dos seguintes dados: 
» U=[1 2 -1; 0 -1 2; 0 0 11]; b=[-3 4 22]’; 
» x = uxb(U,b) 
x = 
 -1 0 2 
Em que ������€ = 
−102 € corresponde à solução do sistema. 
 
1 – Sistemas de Equações não lineares – Métodos iterativos 
 
Considere-se o seguinte sistema de equações lineares Ax=b: 
 
39 
 
„ 3 0 1 −1−1 −5 0 12 1 7 −1−1 0 2 4 … „
�������"… = „
2−72121… 
Cuja matriz do sistema é de diagonal estritamente dominante por linhas, o que leva a 
que neste caso os métodos iterativos de Jacobi e Gauss-Seidel sejam convergentes 
para a solução, para qualquer estimativa inicial do vetor ��8�. 
Função jacobi 
x=jacobi(A,x0,b,maxiter,tol) 
 
Dada um sistema de equações lineares Ax=b, jacobi calcula um valor aproximado do 
vetor solução x, pelo método iterativo de Jacobi. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» x = jacobi(A,x0,b,maxiter,tol 
 
 
 
A – Matriz do sistema, 
x0 – Vetor estimativa inicial, 
b – Vetor dos termos independentes, 
maxiter – Número máximo de iterações, 
tol – Critério de paragem. O processo iterativo termina quando uma determinada norma 
da diferença entre duas iterações consecutivas seja inferior ao valor especificado. A 
norma que está implementada é a norma da soma. 
 
x – valor aproximado da solução obtido na última iteração. 
 
Exemplo: 
Para o sistema dado, com uma estimativa inicial de ��8� = [0 0 0 0]†, 
considerando um número máximo de iterações 100 e como critério de paragem para a 
norma de soma entre duas iterações consecutivas o valor de 0.001, a solução 
aproximada é obtida com a introdução dos seguintes dados: 
» x0=[0 0 0 0]'; maxiter=100; tol=0.001;x=jacobi(A,x0,b,maxiter,tol) 
 (9 iterações) 
x = 
 1.0001 
 2.0000 
 3.0001 
 4.0002 
40 
 
Em que 
‡ˆˆ
ˆˆ‰
���0����0����0��"�0�Š‹‹
‹‹Œ = „1.00012.00003.00014.0002… é a solução aproximada na 9
a
 iteração. 
Se quiséssemos obter uma aproximação melhor, teríamos de alterar o valor de tol, por 
exemplo para um valor de 10-12 , o que representa uma precisão bastante próxima da 
precisão máxima. Assim: 
» x0=[0 0 0 0]'; maxiter=100; tol=1e-12;x=jacobi(A,x0,b,maxiter,tol) 
 (30 iterações) 
x = 
 
 1.0000 
 2.0000 
 3.0000 
 4.0000 
Em que que 
‡ˆˆ
ˆˆ‰
����8�����8�����8��"��8�Š‹‹
‹‹Œ = „1.00002.00003.00004.0000… é a solução aproximada ao fim de 30 iterações. 
De fato, o valor exato para a solução do sistema é: � = [1 2 3 4]† 
 
Função gauss_seidel 
x=gauss_seidel(A,x0,b,maxiter,tol) 
 
Dada um sistema de equações lineares Ax=b, gauss_seidel calcula um valor aproximado 
do vetor solução x, pelo método iterativo de Gauss-Seidel. 
 
A chamada à função, em ambiente Matlab ou Octave, é feita da seguinte forma: 
 
» x = gauss_seidel(A,x0,b,maxiter,tol 
 
 
 
A – Matriz do sistema, 
x0 – Vetor estimativa inicial, 
b – Vetor dos termos independentes, 
maxiter – Número máximo de iterações, 
tol – Critério de paragem. O processo iterativo termina quando uma determinada norma 
da diferença entre duas iterações consecutivas seja inferior ao valor especificado. A 
norma que está implementada é a norma da soma. 
41 
 
 
x – valor aproximado da solução obtido na última iteração. 
 
Exemplo: 
Para o sistema dado, com uma estimativa inicial de ��8� = [0 0 0 0]†, 
considerando um número máximo de iterações 100 e como critério de paragem para a 
norma de soma entre duas iterações consecutivas o valor de 0.001, a solução 
aproximada é obtida com a introdução dos seguintes dados: 
» x0=[0 0 0 0]'; maxiter=100; tol=0.001;x=gauss_seidel(A,x0,b,maxiter,tol) 
 (6 iterações) 
x = 
 1.0000 
 2.0000 
 3.0000 
 4.0000 
Em que 
‡ˆˆ
ˆˆ‰
���$����$����$��"�$�Š‹‹
‹‹Œ = „1.00002.00003.00004.0000… é a solução aproximada na 6
a
 iteração. 
O que reduz o número de iterações relativamente ao método de Jacobi. 
Se quiséssemos obter uma aproximação melhor, em semelhança ao eu fizemos para o 
método anterior de Jacobi, teríamos de alterar o valor de tol, por exemplo para um valor 
de 10-12 , o que representa uma precisão bastante próxima da precisão máxima. Assim: 
» x0=[0 0 0 0]'; maxiter=100; tol=1e-12;x=g_seidel(A,x0,b,maxiter,tol) 
 (15 iterações) 
x = 
 1.0000 
 2.0000 
 3.0000 
 4.0000 
Em que 
‡ˆˆ
ˆˆ‰
����������������"���Š‹‹
‹‹Œ = „1.00002.00003.00004.0000… é a solução aproximada ao fim de 15 iterações. Uma 
redução significativa no número de iterações face ao método de Jacobi. 
 
 
42 
 
2 – Sistemas de Equações Não Lineares 
Considere-se o sistema de equações não lineares: 
f(x)=0 
em que 0 é um vetor nulo, cujo vetor solução aproximada do vetor � = [�� �� … ��]† se pretende aproximar, utilizando o método de Newton. 
O sistema de equações não lineares f(x)=0 de n funções não lineares a n incógnitas é 
representado pelo sistema de equações não lineares: 
Ž!����, ��…��� = 0	!����, ��…��� = 0⋮!����, ��…��� = 0 
A solução aproximada consiste em resolver de uma forma iterativa os seguintes 
sistemas de equações lineares: 
( ) ( )



∆+=
−=∆
+ )()()1(
)()()(
kkk
kkk
xxx
xfxxJ
 
com uma estimativa inicial o vetor ��8�, terminando quando ‘��&'�� − ��&�‘’ ≤ “. 
J(x) é a matriz Jacobiana do vetor de funções f(x), dada por: 
 
J(x)		=
‡ˆˆ
ˆˆˆ
ˆ‰”!� ”��• ”!� ”��• … ”!� ”��•”!� ”��• ”!� ”��• … ”!� ”��•⋮ ⋮ ⋱ ⋮”!� ”��• ”!� ”��• … ”!� ”��• Š‹‹
‹‹‹
‹Œ
 
 
Manual das funções 
 Soleqnl 
 Soleqnl_inter 
43 
 
Função soleqnl 
x=soleqnl(x0,maxiter,tol) 
Dado o sistema de equações não lineares f(x)=0, soleqnl determina uma solução 
aproximada pelo método iterativo de Newton, para uma estimativa inicial x0, e termina 
quando é atingido um número máximo de iterações maxiter, ou quando uma determinada 
norma entre duas interacções consecutivas seja inferior a um determinado valor tol. A 
norma implementada é a norma de soma. 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» x = soleqnl(x0,maxiter,tol) 
 
x0 – vetor com estimativas inicias das variáveis, 
maxiter – número máximo de iterações, 
tol – critério de paragem, como sendo a norma de soma entre duas iterações consecutivas 
 
x – vetor solução aproximado. 
 
A implementação do método de Newton para sistemas de equações não lineares, 
encontra-se no ficheiro soleqnl.m. 
 
Exemplo: 
considere-se o sistema de equações não lineares: 
—������ + �� + �� = −11�� + ������ + �� = −10�� + �� + ������ = −9 
Cuja solução aproximada se pretende obter pelo método de Newton para sistemas de 
equações não lineares. 
O vetor f(x) é dado por: 
 f(x)= ������ + �� + �� + 11�� + ������ + �� + 10�� + �� + ������ + 9 € 
e a sua Jacobiana J(x): 
 
44 
 
J(x)=  ���� ���� + 1 ���� + 1���� + 1 ���� ���� + 1���� + 1 ���� + 1 ���� € 
 
Primeiro é necessário definir a implementação da função f(x) dada, e da sua Jacobiana 
J(x), acedendo ao ficheiro soleqnl.m, e no corpo da função soleqnl já existente no 
ficheiro, introduzir as linhas logo após o comentário: 
 % Introduzir aqui a função e a sua Jacobiana 
 
 f(1) = x(1)*x(2)*x(3) + x(2) + x(3) - 11; 
 f(2) = x(1) + x(1)*x(2)*x(3) + x(3) - 10;f(3) = x(1) + x(2) + x(1)*x(2)*x(3) - 9; 
 J(1,1)=x(2)*x(3); J(1,2)=x(1)*x(3) + 1; J(1,3)=x(1)*x(2) + 1; 
 J(2,1)=x(2)*x(3)+1; J(2,2)=x(1)*x(3); J(2,3)=x(1)*x(2) + 1; 
 J(3,1)=x(2)*x(3)+1; J(3,2)=x(1)*x(3) + 1; J(3,3)=x(1)*x(2); 
 
 
tendo o cuidado de comentar as funções eventualmente existentes. 
Para o sistema dado, com uma estimativa inicial de ��8� = [0 0 0]†, considerando um 
número máximo de iterações 100 e como critério de paragem para a norma de soma entre 
duas iterações consecutivas o valor de 0.004, a solução aproximada é obtida com a 
introdução dos seguintes dados: 
» x0=[0 0 0]'; maxiter=100; tol=0.001;x=soleqnl(x0,maxiter,tol) 
 (7 iterações) 
x = 
 1.0000 
 2.0000 
 3.0000 
Em que „���˜����˜����˜�… = 
1.00002.00003.0000€ é a solução aproximada na 7a iteração. 
 
Função soleqnl_inter 
x=soleqnl_inter(x) 
 
Dado o sistema de equações não lineares f(x)=0, soleqnl realiza uma iteração pelo 
método iterativo de Newton, para uma estimativa x, apresentando, para essa iteração, 
45 
 
também os valores numéricos do vetor função f(x), da matriz Jacobiana J(x), do vetor 
solução ∆x e ainda da norma de soma deste vetor ∆x. 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» x1 = soleqnl_inter(x0) 
 
x0 – vetor com uma estimativa da solução para as variáveis, 
 
x1 – vetor solução aproximado, obtido com uma iteração. 
 
A implementação segue os mesmos passos que a função soleqnl, descrita atrás. 
 
Assim, para o mesmo exemplo da função soleqnl, vão sendo realizadas 6 iterações de 
uma forma interativa e onde podem ser observados os valores numéricos do vetor função 
f(x), da matriz Jacobiana J(x), do vetor solução ∆x e ainda da norma de soma deste vetor 
∆x. De reparar que na primeira chamada da função soleqnl é necessário definir o vetor 
estimativa inicial. Para as restantes iterações o valor do vetor de entrada é o valor 
calculado na iteração anterior, existente em memória. 
» x=[0 0 0]'; [J,f,dx,x,norm]=soleqnl_inter(x) 
J = 
 0 1 1 
 1 0 1 
 1 1 0 
f = 
 -11 -10 -9 
dx = 
 4 
 5 
 6 
x = 
 4 
 5 
 6 
norm = 
 6 
» [J,f,dx,x,norm]= soleqnl_inter (x) 
J = 
 30 25 21 
 31 24 21 
 31 25 20 
46 
 
f = 
 120 120 120 
dx = 
 -1.5789 
 -1.5789 
 -1.5789 
x = 
 2.4211 
 3.4211 
 4.4211 
norm = 
 1.5789 
» [J,f,dx,x,norm]= soleqnl_inter (x) 
J = 
 15.1247 11.7036 9.2825 
 16.1247 10.7036 9.2825 
 16.1247 11.7036 8.2825 
f = 
 33.4597 33.4597 33.4597 
dx = 
 
 -0.9266 
 -0.9266 
 -0.9266 
x = 
 1.4945 
 2.4945 
 3.4945 
norm = 
 0.9266 
» [J,f,dx,x,norm]= soleqnl_inter (x) 
J = 
 8.7168 6.2224 4.7279 
 9.7168 5.2224 4.7279 
 9.7168 6.2224 3.7279 
f = 
 8.0160 8.0160 8.0160 
dx = 
 -0.4076 
 -0.4076 
47 
 
 -0.4076 
x = 
 1.0869 
 2.0869 
 3.0869 
norm = 
 0.4076 
» [J,f,dx,x,norm]= soleqnl_inter (x) 
J = 
 6.4420 4.3551 3.2682 
 7.4420 3.3551 3.2682 
 7.4420 4.3551 2.2682 
f = 
 1.1755 1.1755 1.1755 
dx = 
 -0.0836 
 -0.0836 
 -0.0836 
x = 
 1.0033 
 2.0033 
 3.0033 
norm = 
 0.0836 
» [J,f,dx,x,norm]= soleqnl_inter (x) 
J = 
 6.0166 4.0133 3.0100 
 7.0166 3.0133 3.0100 
 7.0166 4.0133 2.0100 
f = 
 0.0431 0.0431 0.0431 
dx = 
 -0.0033 
 -0.0033 
 -0.0033 
x = 
 1.0000 
 2.0000 
 3.0000 
norm = 
48 
 
 0.0033 
Em que „���$����$����$�… = 
1.00002.00003.0000€ é a solução aproximada na 6a iteração. Pode ver-se ao longo 
das iterações que o valor da norma foi diminuindo, começando no valor 6 e terminando 
no valor 0.0033, o que é demonstrativo que o processo é convergente para a raiz. 
 
49 
 
 
Capítulo VII – Solução de equações diferenciais ordinárias de 
valor inicial 
 
 
Manual das Funções: 
 
 Soleqdif 
 
 
Função Soleqdif 
[xs,ys] = soleqdif(a,b,y0,h,op) 
 
Dada a condição inicial f(a)=y0 da função f(x) que se pretende determinar como solução 
aproximada da equação diferencial y’=f(x,y), soleqdif determina essa solução aproximada 
de f(x), no intervalo ]a, b] com passo h, por diferentes métodos, consoante a opção op 
especificada: 
A chamada à função é feita, em ambiente Matlab ou Octave, digitando na linha de 
comandos: 
 
» [xs, ys] = soleqdif(a,b,y0,h,op) 
 
a – extremo esquerdo do intervalo onde se pretende determinar a solução aproximada da 
equação diferencial, 
b – extremo direito do intervalo onde se pretende determinar a solução aproximada da 
equação diferencial. 
y0 – valor inicial da ordenada f(x) no extremo esquerdo do intervalo, como condição 
inicial, 
h – passo de integração, 
op – opção: 1 – Euler, 2 - Taylor de ordem 2, 3 - Taylor de ordem 3, 4 - Runge Kutta de 
ordem 2 e 5 - Runge Kutta de ordem 3 
 
xs – vetor dos valores das abcissas onde se pretende calcular o valor aproximado da 
solução da equação diferencial, 
ys – vetor dos valores aproximados da solução da função solução da equação diferencial 
nos pontos dados em xs, sendo o primeiro valor, o valor da condição inicial. 
 
A implementação dos métodos de Euler, Taylor de ordem 2, Taylor de ordem 3, Runge 
Kutta de ordem 2 e Runge Kutta de ordem 4, encontra-se nos ficheiros euler.m, 
taylor2.m, taylor3.m, rk2.m e rk4.m respectivamente. 
 
 
50 
 
Exemplo: 
Dada a equação diferencial: 
™´ = !��, ™� = 1 − ™/� com a condição inicial ™�2� = 2. 
cuja solução aproximada no intervalo ]2, 3], com passo h=0.1, se pretende determinar, 
utilizando os métodos aprendidos: Euler, Taylor de ordem 2 e 3 e Runge Kutta de ordem 
2 e 4. 
Primeiro é necessário definir a função !��, ™� dada, acedendo aos ficheiros euler.m, 
taylor2.m, taylor3.m, rk2.m e rk4.m e no corpo da função f(x,y) já existentes nos 
ficheiros, introduzir a linha: 
 y = 1 – y/x; 
 
tal que o resultado seja o seguinte: 
 
 
function y = f(x,y) 
 
 y = 1 – y/x; 
 
end 
tendo o cuidado de comentar as funções eventualmente existentes. 
Para os métodos de Taylor de ordem 2 e 3, terão ainda de ser programadas as segundas e 
terceiras derivadas, respectivamente, de 	™´ = !��, ™�, ou seja: 
™´´ = !-��, ™� + !œ��, ™� × !��, ™� 
e 
™´´´ = !--��, ™� + 2 × !-œ��, ™� + !œ��, ™�ž × !��, ™� + !œœ��, ™� × !��, ™��+ !-��, ™�!œ��, ™� 
Onde neste caso se terá: 
 !-��, ™� = ™/��	 
 !--��, ™� = −2™/�� 
 !œ��, ™� = −1/� 
 !œœ��, ™� = 0 
 !-œ��, ™� = 1/�� 
51 
 
É agora necessário proceder à programação destas derivadas nos ficheiros taylor2.m e 
taylor3.m e proceder à sua gravação em disco. Assim deve aceder-se ao ficheiro 
taylor2.m e no corpo das funções dfx(x,y) e dfy(x,y), já existentes no ficheiro e introduzir 
as linhas referentes às derivadas em ordem a � e a ™ respectivamente: 
 y = y/x^2; 
 e 
 y = -1/x; 
tal que o resultado seja o seguinte: 
 
 
function y = dfx(x,y) 
 
 y = y/x^2; 
 
end 
e 
function y = dfy(x,y) 
 
 y = -1/x; 
 
end 
Deve agora aceder-se ao ficheiro taylor3.m e no corpo das funções dfx(x,y), dfxx(x,y), 
dfy(x,y), dfyy(x,y) e dfxy(x,y), já existentes no ficheiro, introduzir as linhas referentes às 
derivadas de primeira e segunda ordem em ordem a �, derivadas de primeira e segunda 
ordem em ordem a ™ e a derivada em ordem a � e a ™, respectivamente: 
 y = y/x^2; 
 y = -2y/x^3; 
 y = -1/x; 
 y = 0; 
 y = 1/x^2; 
 
tal que o resultadoseja o seguinte: 
 
 
function y = dfx(x,y) 
 
 y = y/x^2; 
 
end 
52 
 
 
function y = dfxx(x,y) 
 
 y = -2y/x^3; 
 
end 
 
function y = dfy(x,y) 
 
 y = -1/x; 
 
end 
 
function y = dfyy(x,y) 
 
 y = 0; 
 
end 
 
function y = dfxy(x,y) 
 
 y = 1/x^2; 
 
end 
Estamos agora em condições de encontrar uma solução aproximada da equação 
diferencial: 
Ÿ™´ = !��, ™� = 1 − ™/�™�2� = 2 
 
no intervalo ]2,3] com passo de integração h=0.1, ou seja nos nós �  = 2 + 0.1h com h = 1…10, onde se irão utilizar os 5 métodos: Euler, Taylor de ordem 2 e 3 e Runge 
Kutta de ordem 2 e 4. 
Assim, para o efeito tem-se: 
a = 2 
b = 3 
y0 = 2 
h = 0.1 
Calculando pelo método de Euler: 
» a = 2; b = 3; y0 = 2; h = 0.1; 
» [x,y1] = soleqdif(a,b,y0,h,1) 
 x = 
 Columns 1 through 9 
53 
 
 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 
2.7000 2.8000 
 Columns 10 through 11 
 2.9000 3.0000 
 
y1 = 
 Columns 1 through 9 
 2.0000 2.0000 2.0048 2.0136 2.0261 2.0417 2.0600 
2.0808 2.1037 
 Columns 10 through 11 
 2.1286 2.1552 
Calculando pelo método de Taylor de 2ª ordem: 
» [x,y2] = soleqdif(a,b,y0,h,2) 
x = 
 Columns 1 through 8 
 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 
2.7000 
 Columns 9 through 11 
 2.8000 2.9000 3.0000 
y2 = 
 Columns 1 through 8 
 2.0000 2.0025 2.0093 2.0198 2.0337 2.0504 2.0697 
2.0912 
 Columns 9 through 11 
 2.1148 2.1401 2.1672 
Calculando pelo método de Taylor de 3ª ordem: 
» [x,y3] = soleqdif(a,b,y0,h,3) 
x = 
 Columns 1 through 7 
 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 
54 
 
 Columns 8 through 11 
 2.7000 2.8000 2.9000 3.0000 
y3 = 
 
 Columns 1 through 7 
 2.0000 2.0024 2.0091 2.0196 2.0333 2.0500 2.0692 
 Columns 8 through 11 
 2.0907 2.1143 2.1396 2.1666 
 
Calculando pelo método de Runge-Kutta de 2ª ordem: 
» [x,y4] = soleqdif(a,b,y0,h,4) 
x = 
 Columns 1 through 7 
 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 
 Columns 8 through 11 
 2.7000 2.8000 2.9000 3.0000 
y4 = 
 Columns 1 through 7 
 2.0000 2.0024 2.0092 2.0197 2.0334 2.0501 2.0694 
 Columns 8 through 11 
 2.0909 2.1144 2.1398 2.1668 
 
E finalmente calculando pelo método de Runge-Kutta de 4ª ordem: 
» [x,y5] = soleqdif(a,b,y0,h,4) 
x = 
 Columns 1 through 7 
 2.0000 2.1000 2.2000 2.3000 2.4000 2.5000 2.6000 
 Columns 8 through 11 
55 
 
 2.7000 2.8000 2.9000 3.0000 
y5 = 
 Columns 1 through 7 
 2.0000 2.0024 2.0091 2.0196 2.0333 2.0500 2.0692 
 Columns 8 through 11 
 2.0907 2.1143 2.1397 2.1667 
Os resultados podem ser agora visualizados no gráfico, com o comando: 
» plot(x,[y1' y2' y3' y4' y5']); grid 
 
 
 
2 2.1 2.2 2.3 2.4 2.5 2.6 2.7 2.8 2.9 3
2
2.02
2.04
2.06
2.08
2.1
2.12
2.14
2.16
2.18