Pseudocódigo Sistemas lineares

Prévia do material em texto

Universidade Federal Rural do Semiárido 
Disciplina: Cálculo Numérico 
Discente: Maria Kauana Mesquita Rodrigues 
 
 
 
 
 
 
 
 
 
Lista – Unidade 3 
 
 
 
 
 
 
 
Pau dos Ferros/RN 
23 de outubro de 2017 
 Soluções 
1.1 
clear 
 
//A contagem de tempo sera iniciada 
tic() 
 
//Estabelecendo quantidade de pontos 
 
n=5; 
 
// Estabelecendo as funções g1, g2 e g3 
 
deff('[y]=g1(x)','y=x^0'); //O programa funciona corretamente se a função 
constante g(x) = c, sempre for definida por g(x) = c*x^0 
 
deff('[y]=g2(x)','y=x'); 
 
//deff('[y]=g3(x)','y=x^2'); 
 
// Estabelecendo tabela de valores (xk, yk) 
 
n=5; 
 
x(1)=0; x(2)=0.25 ; x(3)=0.50; x(4)=0.75; x(5)=1; //x(6)=2.5; x(7)=3; 
 
y(1)=1; y(2)=1.2840 ; y(3)=1.6487; y(4)=2.1170; y(5)=2.7183; //y(6)=1.9; 
y(7)=1.2; 
 
// Projetando os valores <gk,gl> (aqui serão denotados apenas por gkgl 
 
g1g1=sum(g1(x)'*g1(x)); 
g1g2=sum(g1(x)'*g2(x)); 
//g1g3=sum(g1(x)'*g3(x)); 
g2g1=g1g2; 
g2g2=sum(g2(x)'*g2(x)); 
//g2g3=sum(g2(x)'*g3(x)); 
//g3g1=g1g3; 
//g3g2=g2g3; 
//g3g3=sum(g3(x)'*g3(x)); 
yg1=sum(y'*g1(x)); 
yg2=sum(y'*g2(x)); 
//yg3=sum(y'*g3(x)); 
 
// Projetando matrizes para que o método de cramer seja aplicado e o sistema 
linear seja resolvido 
 
A(1,1)=g1g1; A(1,2)=g1g2; //A(1,3)=g1g3; 
A(2,1)=g2g1; A(2,2)=g2g2; //A(2,3)=g2g3; 
//A(3,1)=g3g1; A(3,2)=g3g2; A(3,3)=g3g3; 
 
A1(1,1)=yg1; A1(1,2)=g1g2; //A1(1,3)=g1g3; 
A1(2,1)=yg2; A1(2,2)=g2g2; //A1(2,3)=g2g3; 
//A1(3,1)=yg3; A1(3,2)=g3g2; A1(3,3)=g3g3; 
 
A2(1,1)=g1g1; A2(1,2)=yg1; //A2(1,3)=g1g3; 
A2(2,1)=g2g1; A2(2,2)=yg2; //A2(2,3)=g2g3; 
//A2(3,1)=g3g1; A2(3,2)=yg3; A2(3,3)=g3g3; 
 
//A3(1,1)=g1g1; A3(1,2)=g1g2; //A3(1,3)=yg1; 
//A3(2,1)=g2g1; A3(2,2)=g2g2;// A3(2,3)=yg2; 
//A3(3,1)=g3g1; A3(3,2)=g3g2; A3(3,3)=yg3; 
 
a(1)=det(A1)/det(A); a(2)=det(A2)/det(A); //a(3)=det(A3)/det(A); 
 
disp('solução do sistema'); 
disp(a); 
 
//disp('Função f(x) obtida pelo M.M.Q.') 
 
deff('[y]=f(x)','y=a(1)+a(2)*x'); 
 
disp('Avaliando erro'); 
E=[f(x)-y]'*[f(x)-y]; 
disp(E); 
 
// plot(x,[y,f(x)]) 
 
plot2d(x,y, -4); 
h=abs(x(n)-x(1))/20; 
X=[x(1):h:x(n)]'; 
plot2d(X,f(X), 2); 
 
// final da contagem de tempo 
toc() 
t=toc() 
 
 
1.2 
// Exemplo de programa para encontrar os coeficientes a1, a2, a2 da função 
f(x)=a1*g1(x)+a2*g2(x)+a3*g3(x) pelo método dos mínimos quadrados, de modo 
que o gráfico de f tenha o melhor ajuste possível ao diagrama de dispersão dos 
pontos (xk, yk) dados. 
 
clear 
 
//Iniciando contagem de tempo 
tic() 
 
//Estabelecendo quantidade de pontos 
 
n=4; 
 
// Estabelecendo as funções g1, g2 e g3 
 
deff('[y]=g1(x)','y=x^0'); //O programa funciona corretamente se a função 
constante g(x) = c, sempre for definida por g(x) = c*x^0 
 
deff('[y]=g2(x)','y=x'); 
 
deff('[y]=g3(x)','y=x^2'); 
 
// Estabelecendo tabela de valores (xk, yk) 
 
n=4; 
 
x(1)=0; x(2)=-1 ; x(3)=1; x(4)=2; //x(5)=50; //x(6)=2.5; x(7)=3; 
 
y(1)=0; y(2)=3 ; y(3)=1; y(4)=5; //y(5)=4; //y(6)=1.9; y(7)=1.2; 
 
// Projetando os valores <gk,gl> (aqui serão denotados apenas por gkgl 
 
g1g1=sum(g1(x)'*g1(x)); 
g1g2=sum(g1(x)'*g2(x)); 
g1g3=sum(g1(x)'*g3(x)); 
g2g1=g1g2; 
g2g2=sum(g2(x)'*g2(x)); 
g2g3=sum(g2(x)'*g3(x)); 
g3g1=g1g3; 
g3g2=g2g3; 
g3g3=sum(g3(x)'*g3(x)); 
yg1=sum(y'*g1(x)); 
yg2=sum(y'*g2(x)); 
yg3=sum(y'*g3(x)); 
 
// Projetando matrizes para aplicar o método de cramer para resolução de 
sistema linear 
 
A(1,1)=g1g1; A(1,2)=g1g2; A(1,3)=g1g3; 
A(2,1)=g2g1; A(2,2)=g2g2; A(2,3)=g2g3; 
A(3,1)=g3g1; A(3,2)=g3g2; A(3,3)=g3g3; 
 
A1(1,1)=yg1; A1(1,2)=g1g2; A1(1,3)=g1g3; 
A1(2,1)=yg2; A1(2,2)=g2g2; A1(2,3)=g2g3; 
A1(3,1)=yg3; A1(3,2)=g3g2; A1(3,3)=g3g3; 
 
A2(1,1)=g1g1; A2(1,2)=yg1; A2(1,3)=g1g3; 
A2(2,1)=g2g1; A2(2,2)=yg2; A2(2,3)=g2g3; 
A2(3,1)=g3g1; A2(3,2)=yg3; A2(3,3)=g3g3; 
 
A3(1,1)=g1g1; A3(1,2)=g1g2; A3(1,3)=yg1; 
A3(2,1)=g2g1; A3(2,2)=g2g2; A3(2,3)=yg2; 
A3(3,1)=g3g1; A3(3,2)=g3g2; A3(3,3)=yg3; 
 
a(1)=det(A1)/det(A); a(2)=det(A2)/det(A); a(3)=det(A3)/det(A); 
 
disp('solução do sistema'); 
disp(a); 
 
//disp('Função f(x) obtida pelo M.M.Q.') 
 
deff('[y]=f(x)','y=a(1)+a(2)*x+a(3)*x^2'); 
 
disp('Avaliando erro'); 
E=[f(x)-y]'*[f(x)-y]; 
disp(E); 
 
// plot(x,[y,f(x)]) 
 
plot2d(x,y, -4); 
h=abs(x(n)-x(1))/20; 
X=[-1:h:x(n)]'; 
plot2d(X,f(X), 2); 
 
// fim da contagem de tempo 
toc() 
t=toc() 
 
 0.15 - 1.05x + 1.75x ² 
 
 
 
2.1 
function [y]=f(x); 
 y=%e^x+(3*x^2)-2; 
endfunction; 
 
r=f(0.5); 
 printf('\...f(%g)=%g',0.5,r); 
 
 r=f(1); 
 printf('\...f(%g)=%g',1,r); 
 
 r=f(1.5); 
 printf('\...f(%g)=%g',1.5,r); 
 
 r=f(2); 
 printf('\...f(%g)=%g',2,r); 
 
 r=f(2.5); 
 printf('\...f(%g)=%g',2.5,r); 
 
 I=(1/(2*0.5))*(-3*f(0.5)+4*f(0.5+0.5)-f(0.5+2*(0.5))); 
 disp(I); 
 
 I=(1/(2*0.5))*(-3*f(1)+4*f(1+0.5)-f(1+2*(0.5))); 
 disp(I); 
 
 I=(1/(2*0.5))*(-3*f(1.5)+4*f(1.5+0.5)-f(1.5+2*(0.5))); 
 disp(I); 
 
 I=(1/(2*0.5))*(-3*f(2)+4*f(2+0.5)-f(2+2*(0.5))); 
 disp(I); 
 
 I=(1/(2*0.5))*(-3*f(2.5)+4*f(2.5+0.5)-f(2.5+2*(0.5))); 
 
 disp(I); 
 
...f(0.5)=0.398721...f(1)=3.71828...f(1.5)=9.23169...f(2)=17.3891...f(2.5)=28.932
5 
 4.4452744 
 8.3828547 
 12.928663 
 18.477271 
 25.679214 
 
 
3.1 
x0=1 f(x0)=2 
x1=1.1 f(x1)=2.972 
x2=1.2 f(x2)=4.096 
x3=1.3 f(x3)=5.384 
 A integral é 0.6245 
A integral é 1.8735 
clear // limpa todas variáveis 
clc // limpa o console 
function [y]=f(x); 
 y=2*x^3+x^2+x-2; 
endfunction; 
a=1 //input("Entre com a: ") 
b=1.3//input("Entre com b: ") 
n=3 //input("Entre com n: ") 
 h=(b-a)/n; 
 xi=a 
 soma=0; 
 soma2=0 
 x=0 
 printf('\nx%g=%g',x,a) 
 c=f(xi) 
 printf(' f(x%g)=%g',x,c) 
 xi=a+h; 
 for i=1:n-1; 
 x=x+1 
 printf('\nx%g=%g',x,xi) 
 c=f(xi) 
 printf(' f(x%g)=%g',x,c) 
 soma=soma+f(xi); 
 xi=xi+h; 
 if(modulo(x,3)~=0) 
 soma2=soma2+f(xi) 
 end 
 end 
 printf('\nx%g=%g',n,b) 
 c=f(xi) 
 printf(' f(x%g)=%g',n,c) 
 
IT=(h/8)*(f(a)+f(b)+2*soma+ 3*soma2); 
printf('\n A integral é %g', IT) 
I2S=(3*h/8)*(f(a)+f(b)+2*soma+ 3*soma2); 
 
printf('\n A integral é %g', I2S) 
 
 
4.1 a 
// programa com metodo de Euler modificado (R-K de primeira ordem) para 
resolução do PVI y'=f(x,y); y(x0)=y0. 
 
clear 
 
//comecando contagem de tempo 
 
tic(); 
 
//Estabelecendo intervalo [a,b] onde será resolvido o PVI 
 
a=0; 
b=1; 
 
// Estabelecendo a função f(x,y) 
 
function Y=f(x, y), Y=x-x*y+3*y; endfunction; 
// Estabelecendo partição do intervalo [a,b] 
 
n=5; 
h=(b-a)/n; 
for i=1:n+1, x(i)=a+h*(i-1); end 
 
// Estabelecendo condição inicial (aqui y(1) é o y0) 
 
y(1)=2; 
 
// Projetando aproximação para y(xi) pelo modo de Euler modificado 
 
for k=1:(n) 
 y(k+1)=y(k)+h*f(x(k),y(k)); 
 disp(f(x(k),y(k))); 
 printf('\n O valor no passo %g é %g', k, f(x(k),y(k))) 
end 
 
//fim contagem de tempo 
toc(); 
t=toc(); 
 
plot2d(x,y,-1); 
 
 
 
 
4.1 b) 
// programa com método de Euler modificado (R-K de segunda ordem) para 
resolução do PVI y'=f(x,y); y(x0)=y0. 
 
clear 
 
//comecando contagem de tempo 
 
tic(); 
 
// Estabelecendo intervalo [a,b] onde será resolvido o PVI 
 
a=0; 
b=1; 
 
// Estabelecendo a função f(x,y) 
 
function Y=f(x, y), Y=x-x*y+3*y, endfunction; 
 
// Estabelecendo partição do intervalo [a,b] 
 
n=5; 
h=(b-a)/n; 
for i=1:n+1, x(i)=a+h*(i-1); end 
 
// Estabelecendo condição inicial (aqui y(1) é o y0) 
 
y(1)=2; 
 
// Projetando aproximaçãopara y(xi) pelo método de Euler modificado 
 
for k=1:(n) 
 K1=f(x(k), y(k)); 
 K2=f(x(k)+h,y(k)+h*K1); 
 y(k+1)=y(k)+h*(K1+K2)/2; 
 printf('\n O valor no passo %g é %g', k, f(x(k),y(k))) 
end 
 
//fim contagem de tempo 
toc(); 
t=toc(); 
 
plot2d(x,y,-1); 
 
 
 
 
 
 
 
 
 
O valor no passo 1 é 6 
 9.16 
 
 O valor no passo 2 é 9.16 
 13.4832 
 
 O valor no passo 3 é 13.4832 
 19.148736 
 
 O valor no passo 4 é 19.1487 
 26.228452 
 O valor no passo 5 é 26.2285