Cálculo Numérico
Implementações de Métodos de Cálculo Numérico em Scilab
Método de Euler de 1ª Ordem
Read this doc on Scribd: Implementação do Método de Euler de 1ª Ordem no Scilab
Implementação do Método de Euler de 1ª Ordem no Scilab rftg.development.googlepages.com Data de criação: 4 de Abril de 2007 Última Actualização: 24 de Abril de 2008 Versão: v0.00 – 21/ABR/2008 Autor: Ricardo Filipe Teixeira Gomes AVISO LEGAL Este documento foi elaborado por Ricardo Filipe Teixeira Gomes, a quem se reservam todos os direitos. © 2008 Ricardo Filipe Teixeira Gomes Este documento encontra-se disponível para consulta e utilização desde que sejam respeitados todos os direitos de autor e/ou propriedade intelectual. A cópia parcial ou integral, através de qualquer tipo de meio, dos textos e imagens disponíveis neste documento encontra-se expressamente proíbida a menos que o utilizador respeite os direitos de autoria e/ou propriedade intelectual, citando para isso convenientemente o documento, e incluindo imperterivelmente uma referência clara à página web do autor: “www.rftg.development.googlepages.com”. O material contido neste documento constitui apenas uma informação de carácter geral baseada em experiências pessoais e não pretende de forma alguma influenciar o leitor sobre qualquer matéria específica. O conteúdo deste documento é fornecido como uma comodidade para os leitores e é constituído apenas por informação não vinculativa. O conteúdo deste documento é fornecido “como está” e não se oferece qualquer garantia sobre o mesmo. O autor do documento declina qualquer responsabilidade em caso de prejuízos que possam ocorrer pelo facto de alguém se basear na informação contida neste documento, uma vez que essa informação é de carácter meramente informativo, não se prometendo ou garantindo que seja precisa, completa e actualizada. O mesmo se aplica ao conteúdo de qualquer referência realizado no mesmo. Quaisquer conflitos decorrentes do uso ou relacionados com este documento, ou respeitantes a direitos de autor e/ou propriedade intelectual sobre materiais que façam parte deste documento deverão ser regidos pela Legislação Portuguesa e sujeitos à jurisdição dos tribunais de Portugal. A leitura deste documento e sua utilização pressupõe a aceitação destas condições. © 2008 Ricardo Filipe Teixeira Gomes. rftg.development.googlepages.com Instituto Politécnico do Porto Instituto Superior de Engenharia do Porto Departamento de Engenharia Electrotécnica Licenciatura em Engenharia Electrotécnica, ramo Automação e Sistemas Disciplina de Matemática Aplicada II Implementação do Método de Euler de 1ª Ordem no Scilab 04-04-2007 Professor/Orientador Prof. José Magalhães Trabalho realizado por: Ricardo Filipe Teixeira Gomes // (c) 2008 Ricardo Gomes - TODOS OS DIREITOS RESERVADOS // http://rftg.development.googlepages.com clear //apaga todas as variáveis definidas g=input('Qual a função (introduza de modo explicitado) ','s'); x0=input('Qual o valor inicial de x ? '); y0=input('Qual o valor inicial de y ? '); h=input('Qual passo de iteração h ? '); xF=input('Qual é o valor final de x pertendido ? '); //Parâmetros usados para teste //g = '2*x*y'; //x0 = 0; //y0 = 1; //xF = 2; //h = 0.2; //FUNÇÃO: function out = f(x,y) out = evstr(g) endfunction //Valores iniciais x(1) = x0; y(1) = y0; //cálculo do comprimento dos vectores (de dados) a apresentar c_vect = (xF-x0) / h; //*************************** Passo h *************************** //cálculo do nº de iterações n_iter = (xF - x0) / h; // que neste caso é igual a 'c_vect' for i = 1 : n_iter x(i+1) = x(i) + h; y(i+1) = y(i) + h * f(x(i),y(i)); end //copia para o vector de saída xsaida = x; y1 = y; //*************************************************************** //*************************** Passo h/2 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/2); for i = 1 : n_iter x(i+1) = x(i) + (h/2); y(i+1) = y(i) + (h/2) * f(x(i),y(i)); end //copia para o vector de saída y2(1) = y0; for i = 1 : c_vect y2(i+1) = y(2*i + 1); end //*************************************************************** //*************************** Passo h/4 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/4); for i = 1 : n_iter x(i+1) = x(i) + (h/4); y(i+1) = y(i) + (h/4) * f(x(i),y(i)); end //copia para o vector de saída y4(1) = y0; for i = 1 : c_vect y4(i+1) = y(4*i + 1); end //*************************************************************** //************************** Cálculo EC1 ************************ EC1(1) = 0; for i = 1 : c_vect EC1(i+1) = y2(i+1) - y1(i+1); end //*************************************************************** //************************** Cálculo EC2 ************************ EC2(1) = 0; for i = 1 : c_vect EC2(i+1) = y4(i+1) - y2(i+1); end //*************************************************************** //*********************** Cálculo da Razão ********************** RAZAO(1) = 0; for i = 1 : c_vect if EC2(i+1) == 0 RAZAO(i+1) = 0; else RAZAO(i+1) = EC1(i+1) / EC2(i+1); end end //*************************************************************** saida=[xsaida y1 y2 y4 EC1 EC2 RAZAO] Executando a aplicação no Scilab: -->scipad(); -->Qual a função (introduza de modo explicitado) -->2*x*y Qual o valor inicial de x ? -->0 Qual o valor inicial de y ? -->1 Qual passo de iteração h ? -->0.2 Qual é o valor final de x pertendido ? -->2 saida = 0. 1. 0.2 1. 0.4 1.08 0.6 1.2528 0.8 1.553472 1. 2.050583 1.2 2.8708163 1.4 4.2488081 1.6 6.6281406 1.8 10.870151 2. 18.696659 1. 1.02 1.124448 1.3358442 1.7056059 2.3346334 3.4179032 5.340132 8.8859797 15.717521 29.498643 1. 1.0302757 1.1482993 1.3824406 1.7951229 2.5106623 3.7769489 6.1035986 10.582204 19.660012 39.092999 0. 0.02 0.044448 0.0830442 0.1521339 0.2840503 0.5470870 1.091324 2.2578391 4.8473703 10.801984 0. 0. 0.0102757 1.94633 0.0238513 1.8635484 0.0465964 1.7822037 0.0895170 1.699497 0.1760290 1.6136569 0.3590457 1.5237253 0.7634665 1.4294326 1.6962242 1.3310971 3.9424907 1.2295198 9.5943562 1.1258686
Implementação do Método de Euler de 1ª Ordem no Scilab rftg.development.googlepages.com Data de criação: 4 de Abril de 2007 Última Actualização: 24 de Abril de 2008 Versão: v0.00 – 21/ABR/2008 Autor: Ricardo Filipe Teixeira Gomes AVISO LEGAL Este documento foi elaborado por Ricardo Filipe Teixeira Gomes, a quem se reservam todos os direitos. © 2008 Ricardo Filipe Teixeira Gomes Este documento encontra-se disponível para consulta e utilização desde que sejam respeitados todos os direitos de autor e/ou propriedade intelectual. A cópia parcial ou integral, através de qualquer tipo de meio, dos textos e imagens disponíveis neste documento encontra-se expressamente proíbida a menos que o utilizador respeite os direitos de autoria e/ou propriedade intelectual, citando para isso convenientemente o documento, e incluindo imperterivelmente uma referência clara à página web do autor: “www.rftg.development.googlepages.com”. O material contido neste documento constitui apenas uma informação de carácter geral baseada em experiências pessoais e não pretende de forma alguma influenciar o leitor sobre qualquer matéria específica. O conteúdo deste documento é fornecido como uma comodidade para os leitores e é constituído apenas por informação não vinculativa. O conteúdo deste documento é fornecido “como está” e não se oferece qualquer garantia sobre o mesmo. O autor do documento declina qualquer responsabilidade em caso de prejuízos que possam ocorrer pelo facto de alguém se basear na informação contida neste documento, uma vez que essa informação é de carácter meramente informativo, não se prometendo ou garantindo que seja precisa, completa e actualizada. O mesmo se aplica ao conteúdo de qualquer referência realizado no mesmo. Quaisquer conflitos decorrentes do uso ou relacionados com este documento, ou respeitantes a direitos de autor e/ou propriedade intelectual sobre materiais que façam parte deste documento deverão ser regidos pela Legislação Portuguesa e sujeitos à jurisdição dos tribunais de Portugal. A leitura deste documento e sua utilização pressupõe a aceitação destas condições. © 2008 Ricardo Filipe Teixeira Gomes. rftg.development.googlepages.com Instituto Politécnico do Porto Instituto Superior de Engenharia do Porto Departamento de Engenharia Electrotécnica Licenciatura em Engenharia Electrotécnica, ramo Automação e Sistemas Disciplina de Matemática Aplicada II Implementação do Método de Euler de 1ª Ordem no Scilab 04-04-2007 Professor/Orientador Prof. José Magalhães Trabalho realizado por: Ricardo Filipe Teixeira Gomes // (c) 2008 Ricardo Gomes - TODOS OS DIREITOS RESERVADOS // http://rftg.development.googlepages.com clear //apaga todas as variáveis definidas g=input('Qual a função (introduza de modo explicitado) ','s'); x0=input('Qual o valor inicial de x ? '); y0=input('Qual o valor inicial de y ? '); h=input('Qual passo de iteração h ? '); xF=input('Qual é o valor final de x pertendido ? '); //Parâmetros usados para teste //g = '2*x*y'; //x0 = 0; //y0 = 1; //xF = 2; //h = 0.2; //FUNÇÃO: function out = f(x,y) out = evstr(g) endfunction //Valores iniciais x(1) = x0; y(1) = y0; //cálculo do comprimento dos vectores (de dados) a apresentar c_vect = (xF-x0) / h; //*************************** Passo h *************************** //cálculo do nº de iterações n_iter = (xF - x0) / h; // que neste caso é igual a 'c_vect' for i = 1 : n_iter x(i+1) = x(i) + h; y(i+1) = y(i) + h * f(x(i),y(i)); end //copia para o vector de saída xsaida = x; y1 = y; //*************************************************************** //*************************** Passo h/2 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/2); for i = 1 : n_iter x(i+1) = x(i) + (h/2); y(i+1) = y(i) + (h/2) * f(x(i),y(i)); end //copia para o vector de saída y2(1) = y0; for i = 1 : c_vect y2(i+1) = y(2*i + 1); end //*************************************************************** //*************************** Passo h/4 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/4); for i = 1 : n_iter x(i+1) = x(i) + (h/4); y(i+1) = y(i) + (h/4) * f(x(i),y(i)); end //copia para o vector de saída y4(1) = y0; for i = 1 : c_vect y4(i+1) = y(4*i + 1); end //*************************************************************** //************************** Cálculo EC1 ************************ EC1(1) = 0; for i = 1 : c_vect EC1(i+1) = y2(i+1) - y1(i+1); end //*************************************************************** //************************** Cálculo EC2 ************************ EC2(1) = 0; for i = 1 : c_vect EC2(i+1) = y4(i+1) - y2(i+1); end //*************************************************************** //*********************** Cálculo da Razão ********************** RAZAO(1) = 0; for i = 1 : c_vect if EC2(i+1) == 0 RAZAO(i+1) = 0; else RAZAO(i+1) = EC1(i+1) / EC2(i+1); end end //*************************************************************** saida=[xsaida y1 y2 y4 EC1 EC2 RAZAO] Executando a aplicação no Scilab: -->scipad(); -->Qual a função (introduza de modo explicitado) -->2*x*y Qual o valor inicial de x ? -->0 Qual o valor inicial de y ? -->1 Qual passo de iteração h ? -->0.2 Qual é o valor final de x pertendido ? -->2 saida = 0. 1. 0.2 1. 0.4 1.08 0.6 1.2528 0.8 1.553472 1. 2.050583 1.2 2.8708163 1.4 4.2488081 1.6 6.6281406 1.8 10.870151 2. 18.696659 1. 1.02 1.124448 1.3358442 1.7056059 2.3346334 3.4179032 5.340132 8.8859797 15.717521 29.498643 1. 1.0302757 1.1482993 1.3824406 1.7951229 2.5106623 3.7769489 6.1035986 10.582204 19.660012 39.092999 0. 0.02 0.044448 0.0830442 0.1521339 0.2840503 0.5470870 1.091324 2.2578391 4.8473703 10.801984 0. 0. 0.0102757 1.94633 0.0238513 1.8635484 0.0465964 1.7822037 0.0895170 1.699497 0.1760290 1.6136569 0.3590457 1.5237253 0.7634665 1.4294326 1.6962242 1.3310971 3.9424907 1.2295198 9.5943562 1.1258686
Descarregar implementação: "EULER_1a.sci"
Método de Taylor de 2ª Ordem
Read this doc on Scribd: Implementação do Método de Taylor de 2ª Ordem no Scilab
Implementação do Método de Taylor de 2ª Ordem no Scilab rftg.development.googlepages.com Data de criação: 9 de Abril de 2007 Última Actualização: 24 de Abril de 2008 Versão: v0.00 – 21/ABR/2008 Autor: Ricardo Filipe Teixeira Gomes AVISO LEGAL Este documento foi elaborado por Ricardo Filipe Teixeira Gomes, a quem se reservam todos os direitos. © 2008 Ricardo Filipe Teixeira Gomes Este documento encontra-se disponível para consulta e utilização desde que sejam respeitados todos os direitos de autor e/ou propriedade intelectual. A cópia parcial ou integral, através de qualquer tipo de meio, dos textos e imagens disponíveis neste documento encontra-se expressamente proíbida a menos que o utilizador respeite os direitos de autoria e/ou propriedade intelectual, citando para isso convenientemente o documento, e incluindo imperterivelmente uma referência clara à página web do autor: “www.rftg.development.googlepages.com”. O material contido neste documento constitui apenas uma informação de carácter geral baseada em experiências pessoais e não pretende de forma alguma influenciar o leitor sobre qualquer matéria específica. O conteúdo deste documento é fornecido como uma comodidade para os leitores e é constituído apenas por informação não vinculativa. O conteúdo deste documento é fornecido “como está” e não se oferece qualquer garantia sobre o mesmo. O autor do documento declina qualquer responsabilidade em caso de prejuízos que possam ocorrer pelo facto de alguém se basear na informação contida neste documento, uma vez que essa informação é de carácter meramente informativo, não se prometendo ou garantindo que seja precisa, completa e actualizada. O mesmo se aplica ao conteúdo de qualquer referência realizado no mesmo. Quaisquer conflitos decorrentes do uso ou relacionados com este documento, ou respeitantes a direitos de autor e/ou propriedade intelectual sobre materiais que façam parte deste documento deverão ser regidos pela Legislação Portuguesa e sujeitos à jurisdição dos tribunais de Portugal. A leitura deste documento e sua utilização pressupõe a aceitação destas condições. © 2008 Ricardo Filipe Teixeira Gomes. rftg.development.googlepages.com Instituto Politécnico do Porto Instituto Superior de Engenharia do Porto Departamento de Engenharia Electrotécnica Licenciatura em Engenharia Electrotécnica, ramo Automação e Sistemas Disciplina de Matemática Aplicada II Implementação do Método de Taylor de 2ª Ordem no Scilab 09-04-2007 Professor/Orientador: Prof. José Magalhães Trabalho realizado por: Ricardo Filipe Teixeira Gomes Devido a não ser utilizada a biblioteca simbólica do Scilab é necessário pedir ao utilizador a 1ª e 2ª derivadas, no entanto este facto não é importante para a demonstração em causa, dado que o intuito é provar a validade, e acréscimo de precisão, do método de Taylor de 2ª ordem relativamente ao método de Euler (de 1ª ordem) na resolução de equações diferencias. A fórmula geral de aplicação do método de Taylor de 2ª ordem implementada foi: Visto que a função possui duas variáveis, para ser obtido y’’ será necessário executar os seguintes cálculos: y' ' = d 2 y dF dx ∂F dy = =× + * 2 dx dx ∂y dx dx Pelo que na verdade será pedido ao utilizador a 1ª derivada da função ∂F dy ∂F = F ( x ) e as 2ªs derivadas parciais em ordem a ‘x’ e em ordem ‘y’ ∂y . ∂x dx De seguida é apresentado o código desenvolvido: // (c) 2008 Ricardo Gomes - TODOS OS DIREITOS RESERVADOS // http://rftg.development.googlepages.com //Implementação do Método de Taylor de 2ª ordem //v0.01 09-04-2007 clear //apaga todas as variáveis definidas g=input('Qual a função real (apenas para efeitos de teóricos)? ','s'); g1=input('Qual 1ª derivada da função (introduza de modo explicitado)? ','s'); g2=input('Qual 2ª derivada parcial da função em x (introduza de modo explicitado)? ','s'); g3=input('Qual 2ª derivada parcial da função em y(introduza de modo explicitado)? ','s'); x0=input('Qual o valor inicial de x ? '); xF=input('Qual é o valor final de x pertendido ? '); y0=input('Qual o valor inicial de y ? '); h=input('Qual passo de iteração h ? '); //Parâmetros usados para teste //g = 'exp(x^2)'; //g1 = '2*x*y'; //g2 = '2*y'; //g3 = '2*x'; //x0 = 0; //y0 = 1; //xF = 2; //h = 0.2; //FUNÇÕES: function out = f(x) out = evstr(g) endfunction function out = f1(x,y) out = evstr(g1) endfunction function out = f2(x,y) out = evstr(g2) endfunction function out = f3(x,y) out = evstr (g3) endfunction //Valores iniciais x(1) = x0; y(1) = y0; //cálculo do comprimento dos vectores (de dados) a apresentar c_vect = (xF-x0) / h; //*************************** Passo h *************************** //cálculo do nº de iterações n_iter = (xF - x0) / h; // que neste caso é igual a 'c_vect' for i = 1 : n_iter x(i+1) = x(i) + h; y(i+1) = y(i) + h * f1(x(i),y(i)) + ((h^2)/2) * ( f2(x(i),y(i)) + f3(x(i),y(i))*f1(x(i),y(i)) ); end //copia para o vector de saída xsaida = x; y1 = y; //*************************************************************** //*************************** Passo h/2 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/2); for i = 1 : n_iter x(i+1) = x(i) + (h/2); y(i+1) = y(i) + (h/2) * f1(x(i),y(i)) + (((h/2)^2)/2) * ( f2(x(i),y(i)) + f3(x(i),y(i))*f1(x(i),y(i)) ); end //copia para o vector de saída y2(1) = y0; for i = 1 : c_vect y2(i+1) = y(2*i + 1); end //*************************************************************** //*************************** Passo h/4 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/4); for i = 1 : n_iter x(i+1) = x(i) + (h/4); y(i+1) = y(i) + (h/4) * f1(x(i),y(i)) + (((h/4)^2)/2) * ( f2(x(i),y(i)) + f3(x(i),y(i))*f1(x(i),y(i)) ); end //copia para o vector de saída y4(1) = y0; for i = 1 : c_vect y4(i+1) = y(4*i + 1); end //*************************************************************** //************************** Cálculo EC1 ************************ EC1(1) = 0; for i = 1 : c_vect EC1(i+1) = y2(i+1) - y1(i+1); end //*************************************************************** //************************** Cálculo EC2 ************************ EC2(1) = 0; for i = 1 : c_vect EC2(i+1) = y4(i+1) - y2(i+1); end //*************************************************************** //*********************** Cálculo da Razão ********************** RAZAO(1) = 0; for i = 1 : c_vect if EC2(i+1) == 0 RAZAO(i+1) = 0; else RAZAO(i+1) = EC1(i+1) / EC2(i+1); end end //*************************************************************** //******************* Cálculo da solução real ************************ for i = 1:c_vect+1 yreal(i) = f(xsaida(i)); end //******************************************************************** printf("\n\n\n\n") printf("*********************************\n") printf("* Tabela método Taylor 2ª ordem *\n") printf("*********************************\n") printf("\nOs parâmetros apresentados são, respectivamente:\n") printf("\n x; y(x); y para h; y para h/2; y para h/4; Erro 1; Erro 2; Razão \n\n") Tabela=[xsaida yreal y1 y2 y4 EC1 EC2 RAZAO] printf("\n\nFim.\n") Executando a aplicação no Scilab: -->scipad() -->Qual a função real (apenas para efeitos de teóricos)? -->exp(x^2) Qual 1ª derivada da função (introduza de modo explicitado)? -->2*x*y Qual 2ª derivada parcial da função em x (introduza de modo explicitado)? -->2*y Qual 2ª derivada parcial da função em y(introduza de modo explicitado)? -->2*x Qual o valor inicial de x ? -->0 Qual é o valor final de x pertendido ? -->2 Qual o valor inicial de y ? -->1 Qual passo de iteração h ? -->0.2 ********************************* * Tabela método Taylor 2ª ordem * ********************************* Os parâmetros apresentados são, respectivamente: x; y(x); y para h; y para h/2; y para h/4; Erro 1; Erro 2; Razão Tabela 0. 0.2 0.4 0.6 0.8 1. 1.2 1.4 1.6 1.8 2. Fim. = 1. 1.0408108 1.1735109 1.4333294 1.8964809 2.7182818 4.2206958 7.0993271 12.935817 25.533722 54.59815 1. 1.04 1.168128 1.4167056 1.8541843 2.6166249 3.9772699 6.5036318 11.42558 21.534934 43.483338 1. 1.040502 1.1718627 1.4284046 1.8839578 2.6878302 4.1464242 6.9133307 12.451365 24.211098 50.794026 1. 1.0407196 1.1730583 1.4319937 1.8930814 2.7099647 4.2002253 7.0474891 12.799052 25.154888 53.490859 0. 0.000502 0.0037347 0.0116989 0.0297735 0.0712052 0.1691543 0.4096990 1.0257847 2.6761647 7.3106883 0. 0.0002176 0.0011956 0.0035891 0.0091236 0.0221345 0.0538011 0.1341584 0.3476866 0.9437900 2.6968328 0. 2.3069436 3.1236791 3.259548 3.2633545 3.2169335 3.1440693 3.053845 2.9503143 2.8355509 2.7108422 Comparando com os valores obtidos através do método de Euler (de 1ª ordem) (ver abaixo) é possível observar que ocorreu um aumento da precisão de cálculo devido à aplicação do método de Taylor de 2ª ordem, visto que houve uma redução significativa do erro do relativamente ao valor analítico. Método de Euler: -->scipad(); -->Qual a função (introduza de modo explicitado) -->2*x*y Qual o valor inicial de x ? -->0 Qual o valor inicial de y ? -->1 Qual passo de iteração h ? -->0.2 Qual é o valor final de x pertendido ? -->2 saida = 0. 0.2 0.4 0.6 0.8 1. 1.2 1.4 1.6 1.8 2. 1. 1. 1.08 1.2528 1.553472 2.050583 2.8708163 4.2488081 6.6281406 10.870151 18.696659 1. 1.02 1.124448 1.3358442 1.7056059 2.3346334 3.4179032 5.340132 8.8859797 15.717521 29.498643 1. 1.0302757 1.1482993 1.3824406 1.7951229 2.5106623 3.7769489 6.1035986 10.582204 19.660012 39.092999 0. 0.02 0.044448 0.0830442 0.1521339 0.2840503 0.5470870 1.091324 2.2578391 4.8473703 10.801984 0. 0.0102757 0.0238513 0.0465964 0.0895170 0.1760290 0.3590457 0.7634665 1.6962242 3.9424907 9.5943562 0. 1.94633 1.8635484 1.7822037 1.699497 1.6136569 1.5237253 1.4294326 1.3310971 1.2295198 1.1258686
Implementação do Método de Taylor de 2ª Ordem no Scilab rftg.development.googlepages.com Data de criação: 9 de Abril de 2007 Última Actualização: 24 de Abril de 2008 Versão: v0.00 – 21/ABR/2008 Autor: Ricardo Filipe Teixeira Gomes AVISO LEGAL Este documento foi elaborado por Ricardo Filipe Teixeira Gomes, a quem se reservam todos os direitos. © 2008 Ricardo Filipe Teixeira Gomes Este documento encontra-se disponível para consulta e utilização desde que sejam respeitados todos os direitos de autor e/ou propriedade intelectual. A cópia parcial ou integral, através de qualquer tipo de meio, dos textos e imagens disponíveis neste documento encontra-se expressamente proíbida a menos que o utilizador respeite os direitos de autoria e/ou propriedade intelectual, citando para isso convenientemente o documento, e incluindo imperterivelmente uma referência clara à página web do autor: “www.rftg.development.googlepages.com”. O material contido neste documento constitui apenas uma informação de carácter geral baseada em experiências pessoais e não pretende de forma alguma influenciar o leitor sobre qualquer matéria específica. O conteúdo deste documento é fornecido como uma comodidade para os leitores e é constituído apenas por informação não vinculativa. O conteúdo deste documento é fornecido “como está” e não se oferece qualquer garantia sobre o mesmo. O autor do documento declina qualquer responsabilidade em caso de prejuízos que possam ocorrer pelo facto de alguém se basear na informação contida neste documento, uma vez que essa informação é de carácter meramente informativo, não se prometendo ou garantindo que seja precisa, completa e actualizada. O mesmo se aplica ao conteúdo de qualquer referência realizado no mesmo. Quaisquer conflitos decorrentes do uso ou relacionados com este documento, ou respeitantes a direitos de autor e/ou propriedade intelectual sobre materiais que façam parte deste documento deverão ser regidos pela Legislação Portuguesa e sujeitos à jurisdição dos tribunais de Portugal. A leitura deste documento e sua utilização pressupõe a aceitação destas condições. © 2008 Ricardo Filipe Teixeira Gomes. rftg.development.googlepages.com Instituto Politécnico do Porto Instituto Superior de Engenharia do Porto Departamento de Engenharia Electrotécnica Licenciatura em Engenharia Electrotécnica, ramo Automação e Sistemas Disciplina de Matemática Aplicada II Implementação do Método de Taylor de 2ª Ordem no Scilab 09-04-2007 Professor/Orientador: Prof. José Magalhães Trabalho realizado por: Ricardo Filipe Teixeira Gomes Devido a não ser utilizada a biblioteca simbólica do Scilab é necessário pedir ao utilizador a 1ª e 2ª derivadas, no entanto este facto não é importante para a demonstração em causa, dado que o intuito é provar a validade, e acréscimo de precisão, do método de Taylor de 2ª ordem relativamente ao método de Euler (de 1ª ordem) na resolução de equações diferencias. A fórmula geral de aplicação do método de Taylor de 2ª ordem implementada foi: Visto que a função possui duas variáveis, para ser obtido y’’ será necessário executar os seguintes cálculos: y' ' = d 2 y dF dx ∂F dy = =× + * 2 dx dx ∂y dx dx Pelo que na verdade será pedido ao utilizador a 1ª derivada da função ∂F dy ∂F = F ( x ) e as 2ªs derivadas parciais em ordem a ‘x’ e em ordem ‘y’ ∂y . ∂x dx De seguida é apresentado o código desenvolvido: // (c) 2008 Ricardo Gomes - TODOS OS DIREITOS RESERVADOS // http://rftg.development.googlepages.com //Implementação do Método de Taylor de 2ª ordem //v0.01 09-04-2007 clear //apaga todas as variáveis definidas g=input('Qual a função real (apenas para efeitos de teóricos)? ','s'); g1=input('Qual 1ª derivada da função (introduza de modo explicitado)? ','s'); g2=input('Qual 2ª derivada parcial da função em x (introduza de modo explicitado)? ','s'); g3=input('Qual 2ª derivada parcial da função em y(introduza de modo explicitado)? ','s'); x0=input('Qual o valor inicial de x ? '); xF=input('Qual é o valor final de x pertendido ? '); y0=input('Qual o valor inicial de y ? '); h=input('Qual passo de iteração h ? '); //Parâmetros usados para teste //g = 'exp(x^2)'; //g1 = '2*x*y'; //g2 = '2*y'; //g3 = '2*x'; //x0 = 0; //y0 = 1; //xF = 2; //h = 0.2; //FUNÇÕES: function out = f(x) out = evstr(g) endfunction function out = f1(x,y) out = evstr(g1) endfunction function out = f2(x,y) out = evstr(g2) endfunction function out = f3(x,y) out = evstr (g3) endfunction //Valores iniciais x(1) = x0; y(1) = y0; //cálculo do comprimento dos vectores (de dados) a apresentar c_vect = (xF-x0) / h; //*************************** Passo h *************************** //cálculo do nº de iterações n_iter = (xF - x0) / h; // que neste caso é igual a 'c_vect' for i = 1 : n_iter x(i+1) = x(i) + h; y(i+1) = y(i) + h * f1(x(i),y(i)) + ((h^2)/2) * ( f2(x(i),y(i)) + f3(x(i),y(i))*f1(x(i),y(i)) ); end //copia para o vector de saída xsaida = x; y1 = y; //*************************************************************** //*************************** Passo h/2 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/2); for i = 1 : n_iter x(i+1) = x(i) + (h/2); y(i+1) = y(i) + (h/2) * f1(x(i),y(i)) + (((h/2)^2)/2) * ( f2(x(i),y(i)) + f3(x(i),y(i))*f1(x(i),y(i)) ); end //copia para o vector de saída y2(1) = y0; for i = 1 : c_vect y2(i+1) = y(2*i + 1); end //*************************************************************** //*************************** Passo h/4 ************************* //cálculo do nº de iterações n_iter = (xF - x0) / (h/4); for i = 1 : n_iter x(i+1) = x(i) + (h/4); y(i+1) = y(i) + (h/4) * f1(x(i),y(i)) + (((h/4)^2)/2) * ( f2(x(i),y(i)) + f3(x(i),y(i))*f1(x(i),y(i)) ); end //copia para o vector de saída y4(1) = y0; for i = 1 : c_vect y4(i+1) = y(4*i + 1); end //*************************************************************** //************************** Cálculo EC1 ************************ EC1(1) = 0; for i = 1 : c_vect EC1(i+1) = y2(i+1) - y1(i+1); end //*************************************************************** //************************** Cálculo EC2 ************************ EC2(1) = 0; for i = 1 : c_vect EC2(i+1) = y4(i+1) - y2(i+1); end //*************************************************************** //*********************** Cálculo da Razão ********************** RAZAO(1) = 0; for i = 1 : c_vect if EC2(i+1) == 0 RAZAO(i+1) = 0; else RAZAO(i+1) = EC1(i+1) / EC2(i+1); end end //*************************************************************** //******************* Cálculo da solução real ************************ for i = 1:c_vect+1 yreal(i) = f(xsaida(i)); end //******************************************************************** printf("\n\n\n\n") printf("*********************************\n") printf("* Tabela método Taylor 2ª ordem *\n") printf("*********************************\n") printf("\nOs parâmetros apresentados são, respectivamente:\n") printf("\n x; y(x); y para h; y para h/2; y para h/4; Erro 1; Erro 2; Razão \n\n") Tabela=[xsaida yreal y1 y2 y4 EC1 EC2 RAZAO] printf("\n\nFim.\n") Executando a aplicação no Scilab: -->scipad() -->Qual a função real (apenas para efeitos de teóricos)? -->exp(x^2) Qual 1ª derivada da função (introduza de modo explicitado)? -->2*x*y Qual 2ª derivada parcial da função em x (introduza de modo explicitado)? -->2*y Qual 2ª derivada parcial da função em y(introduza de modo explicitado)? -->2*x Qual o valor inicial de x ? -->0 Qual é o valor final de x pertendido ? -->2 Qual o valor inicial de y ? -->1 Qual passo de iteração h ? -->0.2 ********************************* * Tabela método Taylor 2ª ordem * ********************************* Os parâmetros apresentados são, respectivamente: x; y(x); y para h; y para h/2; y para h/4; Erro 1; Erro 2; Razão Tabela 0. 0.2 0.4 0.6 0.8 1. 1.2 1.4 1.6 1.8 2. Fim. = 1. 1.0408108 1.1735109 1.4333294 1.8964809 2.7182818 4.2206958 7.0993271 12.935817 25.533722 54.59815 1. 1.04 1.168128 1.4167056 1.8541843 2.6166249 3.9772699 6.5036318 11.42558 21.534934 43.483338 1. 1.040502 1.1718627 1.4284046 1.8839578 2.6878302 4.1464242 6.9133307 12.451365 24.211098 50.794026 1. 1.0407196 1.1730583 1.4319937 1.8930814 2.7099647 4.2002253 7.0474891 12.799052 25.154888 53.490859 0. 0.000502 0.0037347 0.0116989 0.0297735 0.0712052 0.1691543 0.4096990 1.0257847 2.6761647 7.3106883 0. 0.0002176 0.0011956 0.0035891 0.0091236 0.0221345 0.0538011 0.1341584 0.3476866 0.9437900 2.6968328 0. 2.3069436 3.1236791 3.259548 3.2633545 3.2169335 3.1440693 3.053845 2.9503143 2.8355509 2.7108422 Comparando com os valores obtidos através do método de Euler (de 1ª ordem) (ver abaixo) é possível observar que ocorreu um aumento da precisão de cálculo devido à aplicação do método de Taylor de 2ª ordem, visto que houve uma redução significativa do erro do relativamente ao valor analítico. Método de Euler: -->scipad(); -->Qual a função (introduza de modo explicitado) -->2*x*y Qual o valor inicial de x ? -->0 Qual o valor inicial de y ? -->1 Qual passo de iteração h ? -->0.2 Qual é o valor final de x pertendido ? -->2 saida = 0. 0.2 0.4 0.6 0.8 1. 1.2 1.4 1.6 1.8 2. 1. 1. 1.08 1.2528 1.553472 2.050583 2.8708163 4.2488081 6.6281406 10.870151 18.696659 1. 1.02 1.124448 1.3358442 1.7056059 2.3346334 3.4179032 5.340132 8.8859797 15.717521 29.498643 1. 1.0302757 1.1482993 1.3824406 1.7951229 2.5106623 3.7769489 6.1035986 10.582204 19.660012 39.092999 0. 0.02 0.044448 0.0830442 0.1521339 0.2840503 0.5470870 1.091324 2.2578391 4.8473703 10.801984 0. 0.0102757 0.0238513 0.0465964 0.0895170 0.1760290 0.3590457 0.7634665 1.6962242 3.9424907 9.5943562 0. 1.94633 1.8635484 1.7822037 1.699497 1.6136569 1.5237253 1.4294326 1.3310971 1.2295198 1.1258686
Descarregar implementação: "TAYLOR_2a.sci"
Método de Runge-Kutta de 2ª e 4ª ordem
Descarregar implementação: "RK2.sci"; "RK4.sci"
Método de Euler de 1ª Ordem para resolução de Equações Diferencias Ordinárias de 2ª Ordem
Descarregar implementação: "EULER_1a_ODEs_2a.sci"
Método de Taylor de 2ª Ordem para resolução de Equações Diferencias Ordinárias de 2ª Ordem
Descarregar implementação: "TAYLOR_2a_ODEs_2a.sci"
Método de Runge-Kutta de 2ª Ordem para resolução de Equações Diferencias Ordinárias de 2ª Ordem
Descarregar implementação: "RK2_ODEs_2a.sci"