Equações

VictorCacciari · Março 20, 2009

Boas pessoal!

hoje, estava dando contínuidade a um pequeno projeto, e me deparei com um problemão.

É um interpretador matemático.

Já tenho quase tudo pronto, funções, gráficos com OpenGL, variáveis, etc...

Mas não consigo de forma alguma resolver equações de funções descontínuas em R (logarítmicas e fracionárias)

(Para as funções contínuas, eu uso o método de Newton e resolve-se fácil.)

O que eu pensei.

Imaginemos uma equação:

x^2 + 13 = log (x/2)

Eu conheço duas possibilidades de encontrar as raízes:

-isolar o 'x'

-separar as funções e encontrar o ponto de intersecção

Eu tenho que ir pela segunda, pois eu trato as expresões e funções como postfix, o que deixa meio complicado de isolar o 'x'

então teríamos duas funções:

f1(x) = x^2+13
f2(x) = log(x/2)

Eu não consigo encontrar o ponto de intersecção dos gráficos...

alguém tem alguma idéia?

alguém conhece outro método?

Obrigado!!

Rarylson Freitas · Maio 6, 2009

No caso do Newton-Rapson, acho que o problema está no teu valor inicial, jah que Newton-Rapson pode divergir e devolver um x intermediario fora dos limites em que as funções são definidas. No meu programa eu nem tratei isso!!!

Agora uma função que com certeza vai lhe dar as raízes é o método da bisserção, de eficiencia ruim, porém que sempre resolve, desde que você inicie suas iterações com dois valores tais que a solução lá se encontre e que o numero de solucoes no intervalo seja impar. (esse método usa 2 valores iniciais, e não 1 como o de Newton-Raphson)

Qq coisa dá uma estuda em um livro de cálculo numérico

Vai aí um programa exemplo:

#include <stdio.h>
#include <math.h>

///////////////////////////////
/* definicoes das funcoes */

/* funcao f(x) - a equacao a ser resolvida eh f(x) = 0. caso tenhamos h(x) = g(x), definir f(x) = h(x)-g(x)  */

float f (float x) {
	  return x - cos(x);	  // para testar outras funcoes, alterar esta e recompilar o programa
}

/* derivada de f(x) */

float f_linha (float x, float dx) {
	  return (f(x+dx)-f(x-dx)) / (2*dx);
}

/* calcula uma raiz pelo metodo de Newton-Raphson */
void newtonraph (void) {
	int nmi, k;
	float x0, x, tol, dx;
	//leitura das variaveis de interesse
	printf ("Digite o valor inicial de x:\n");
	scanf ("%f", &x0);	 
	printf ("Digite a tolerancia requerida:\n");
	scanf ("%f", &tol);
	printf ("Digite o numero maximo de interacoes:\n");
	scanf ("%d", &nmi);
	//calculo da raiz
	k=0;
	dx=tol;
	x = x0 - f(x0)/f_linha(x0, dx);
	printf ("k\tx\tf(x)\tf'(x)\n");   // cabecalho da tabela
	printf ("%d\t%.4f\t%.4f\t%.4f\n", 0, x0, f(x0), f_linha(x0, dx)); // tabela - 1a linha
	printf ("%d\t%.4f\t%.4f\t%.4f\n", 1, x, f(x), f_linha(x, dx)); // tabela - 2a linha	
	while( ( fabs(x-x0)>tol || fabs(f(x))>tol ) && k<=nmi ) {
		 x0 = x;
		 x = x0 - f(x0)/f_linha(x0, dx);
		 printf ("%d\t%.4f\t%.4f\t%.4f\n", k+2, x, f(x), f_linha(x, dx)); // tabela
		 k++;
	}
	if (k>nmi)
		printf("Ultrapassou o numero maximo de interacoes\n");
	else
		printf("Raiz: %f\n", x);
}

/* calcula uma raiz pelo metodo da bissecao */
void bissecao (void) {
	int nmi, k;
	float xi, xf, xm, tol;
	//leitura das variaveis de interesse
	printf ("Digite o valor inicial de x:\n");
	scanf ("%f", ξ);
	printf ("Digite o valor final de x:\n");
	scanf ("%f", &xf);
	printf ("Digite a tolerancia requerida:\n");
	scanf ("%f", &tol);
	printf ("Digite o numero maximo de interacoes:\n");
	scanf ("%d", &nmi);
	 //calculo da raiz
	if ( (f(xi))*(f(xf))<=0 ) {
	  k=0;
	  xm = (xi+xf)/2;
	  printf ("k\txi\tf(xi)\tx(f)\tf(xf)\txm\tf(xm)\n");   // cabecalho da tabela
	  while( ( fabs(xf-xi)>tol || fabs(f(xm))>tol ) && k<=nmi ) {
			 xm = (xi+xf)/2;
			 printf ("%d\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\t%.4f\n", k, xi, f(xi), xf, f(xf), xm, f(xm)); // tabela
			 if ( f(xm)*f(xi)<=0 )
				xf = xm;
			 else
				xi = xm;
			 k++;
	  }
	  if (k>nmi)
		printf("Ultrapassou o numero maximo de interacoes\n");
	  else
		printf("Raiz: %f\n", xm);
	}
	else
		printf("Nao ha garantia de existencia de raiz neste intevalo\n");
}

Esse código eh soh pra estudo (tanto que ele tem um monte de printf pra facilitar o entendimento), na prática você deve fazer um melhor.

VictorCacciari · Maio 6, 2009

Rarylson Freitar

Muito obrigado pelo exemplo.

Vou estudá-lo com calma quando chegar em casa.

Eu não conhecia esse método da bisseção.

Eu estava encontrando as raizes de funções continuas da seguinte forma:

->Pelo corolário do teorema de Bolzano, eu encontrava [a, b] tal que f(a) * f(b) < 0
->Começava as iterações do método de Newton, (normalmente são necessárias 4 iterações para encontrar um x com uma taxa de erro aceitável
->Continuava pela busca de outro intervalo a partir de b, A função poderia ser quadrática por exemplo, então teria mais de uma raiz.
ps.: Se a função fosse composta de pelo menos uma função periódica, eu não buscava mais intervalos.. =P

E o meu problema é quando a funão é descontínua no intervalo [a, b]. Pelo que li do método da bissecção, ele também não me resolve o problema... =/

Obrigado novamente!

abraços

viciado · Maio 12, 2009

VictorCacciari depois de pronto se você puder disponibilizar para o pessoal vai ser uma boa.

VictorCacciari · Maio 13, 2009

Viciado.

Claro que disponibilizarei.

Tou recodificando tudo, usando shared-librarys e adicionando recursão de comandos, espero que dentro de 1 semana esteja completo!

Viva o open-source! =D

Rarylson Freitas · Maio 23, 2009

Tem razão, o método da bisseção usa a hipótese da função ser contínua para achar a raíz (ele usa o Teorema de Bolzano).

O de newton também pode dar erro em situações como essa.

Problemas de cálculo numérico algumas vezes se tornam bem difíceis. Talvez seja o caso de buscar (não estou dizendo que não está fazendo isso) um algoritmo que resolva somente seu caso específico, porque algoritmos gerais são complicados.

Para se ter uma idéia, até o Matlab às vezes não consegue resolver certas equações ou resolve pequenas imprecisões.

Mas o importante é que o geral do cálculo numérico você entendeu.

Se o teu problema for mais complicado, você pode pedir um auxílio a alguém que estuda ou é formado Ciências de Computação e é bom em Análise Numérica.

O curso que eu tive disso foi muito superficial, só mesmo pra ter uma idéia.

VictorCacciari · Maio 23, 2009

Raryson

Obrigado por toda a ajuda.

Eu acabei optando pelo método de Newton, ele converge mais rapidamente no resultado. Fiz uma trava do tipo: "Se ocorrerem mais de 40 iterações, retornar: não encontrado"

Pois em algumas funções (descontinuas e de grau acima de 3) o método de Newton pode entrar num loop infinito.

Vou fazer uma ou duas funções específicas para resolver equações de grau 2 e 3, assim posso meter complexos no meio também.

E na documentação da minha função "solve" eu escrevo que deve ser usada com cuidado, e requer uma pré analise da função por parte do usuário.

O problema é que não conheço ninguém que possa me ajudar no assunto. Ainda não estou na faculdade.

No mais, o projeto ja está quase pronto, falta só administrar multiplas bibliotecas dinamicamente e ja posso soltar a primeira release.

Rarylson Freitas · Maio 27, 2009

Taki um que eu tinha implementado que calcula as raizes de qq polinomio de uma variável.

Mais uma vez, como eu fiz o programa a 1 ano atrás e fiz como um trabalho q era pra ser entregue, ele nao esta modularizado (está solto dentro da main) e tem partes dele que nem eu consigo mais entender, rsrsrs.

O método usado eh o de newton-bairston.

Apesar de naum tah mto bom, o programa da pra ser usado para entender melhor o método ou para comparar os resultados de algum outro q você venha a implementar.

Vale ao menos pra compilar e se divertir, rsrs.

#include <stdio.h>
#include <math.h>
#define TAM 100

/* programa modificado do algoritmo do livro */
/* imprime passo-a-passo */
/* erro consertado: pode ser que p0 e q0 nao convirjam nunca para p e q, isso ocorre se d=0 para p0 e q0;
erro consertado: pode haver divisao por zero no calculo de p0 e q0
solucao: dar oportunidade do usuario alterar manualmente os valores de p0 e q0 */

// calcula as raizes da eq de 2 grau x^2 + p*x + q = 0
void eq2grau (float p, float q, float pr[TAM], float pi[TAM], int *j) {
	 float delta, r, im;
	 
	 delta = p*p - 4*q;
	 r = -p/2;
	 im = sqrt( fabs(delta) ) / 2;
	 if (delta>=0){
		pr[*j] = r + im;
		pr[*j+1] = r - im;
		pi[*j] = pi[*j+1] = 0;   
	 }
	 else {
		  pr[*j] = pr[*j+1] = r;
		  pi[*j] = im;
		  pi[*j+1] = -im;
	 }
	 *j=*j+2;
}

// imprimir eq de grau n
void imprimir_eq (float a[TAM], int n) {
	 int i;
	 
	 if (a[n]==1)
		  printf ("x^%d ", n);
	 else
		  printf ("%f x^%d ", a[n], n);
		 
	 for (i=n-1; i>=0; i--) {
		 if (a[i]>0)
			 printf ("+ %f x^%d ", a[i], i);
		 else
			 printf ("- %f x^%d ", fabs(a[i]), i);;
	 }
	 printf ("= 0\n");
}



int main(void) { 
	float a[TAM], v[TAM], b[TAM+2], c[TAM+2];	// vetores: a => coef; v => coef ordenados; b => recursao e coef da eq fatorada; c => recursao
	float c1, d, t, troca;					   // aux: c1 => corrigir erro; t => transformar a[n]=1
	float p, q, p0, q0, dp, dq;				  // coef da eq de 2 grau
	float tol;
	int n, m, i, j, k, nmi;					  // n => grau da eq; m => grau da eq fatorada; k => cont de interacoes; j => sol[j]
	float pr[TAM], pi[TAM];					  // pr, pi => partes real e imaginaria da solucao
	
	// ler
	printf ("Digite a ordem da equacao:\n");
	scanf ("%d", &n);
	for (i=n; i>=0; i--) {
		printf ("Digite o coeficiente de x^%d:\n", i);
		scanf ("%f", &a[i]);
	}
	printf ("Digite a tolerancia:\n");
	scanf ("%f", &tol);
	printf ("Digite o numero maximo de iteracoes:\n");
	scanf ("%d", &nmi);
	
	printf ("Equacao inicial (grau %d):\n", n);					  // imprimir eq inicial
	imprimir_eq (a, n);
	
	// transformar a[n]=1
	if (a[n]!=1) {
		t=a[n];
		for (i=n; i>=0; i--) {
			a[i] = a[i]/t;
			v[i] = a[i];
		}
		
		printf ("Equacao dividida pelo coef. de x^%d\n", n);		 // imprimir eq transformada
		imprimir_eq (a, n);
	}
	else
		for (i=n; i>=0; i--)
			v[i] = a[i];
	
	// ordenar v[n] para calc p0 e q0
	m=n;
	if (m>2) {
		for (i=0; i<=m-1; i++) {
			for (j=0; j <= m-i-1; j++) {
				if ( fabs(v[j]) > fabs(v[j+1]) ) {
				   troca = v[j];
				   v[j] = v[j+1];
				   v[j+1] = troca;
				}
			}
		}
		if (v[m-2]!=0) {
		   p0 = v[m]/v[m-2];
		   q0 = v[m-1]/v[m-2];
		}
		else {								  // erro: nao calcula p0 e q0; divisao por zero
			 printf ("Nao foi possivel calcular os valores iniciais p0 e q0: divisao por zero\n");
			 printf ("Digite um valor para p0:\n");
			 scanf ("%f", &p0);
			 printf ("Digite um valor para q0:\n");
			 scanf ("%f", &q0);
		}
	}
	j=0;
	k=0;
	// correcoes para grau = 1 ou 2
	b[3] = a[1];
	b[2] = a[0];
	
	// enquanto grau > 2, calcular 2 raizes e fatorar (m = m-2)
	while (m>2 && k<=nmi) {
		 
		 printf ("p0 = %f, q0 = %f\n", p0, q0);		  // imprimir p0 e q0
		 // valores iniciais
		 b[m+2] = b[m+1] = c[m+2] = c[m+1] = 0;
		 p=p0; q=q0;
		 dp = dq = 1;
		 k=0;
		 
		 // calculo de p e q
		 while ( (fabs(dp) + fabs(dq) > tol) && k<=nmi ) {
			   
			   for (i=m; i>=1; i--) {
				   b[i] = a[i] - p*b[i+1] - q*b[i+2];
				   c[i] = b[i] - p*c[i+1] - q*c[i+2];
			   }
			   
			   b[0] = a[0] - p*b[1] - q*b[2];
			   c1 = c[1] - b[1];
			   d = c[2]*c[2] - c1*c[3];
			   
			   if (fabs(d)>0) {
				   dp = (b[1]*c[2] - b[0]*c[3]) / d;
				   dq = (b[0]*c[2] - b[1]*c1) / d;
				   p = p+dp;
				   q = q+dq;
			   }
			   // erro: nao convergiu com p0 e q0
			   else {
				   printf ("Nao convergiu com os valores iniciais p0 e q0\n");
				   printf ("Digite um novo valor para p0:\n");
				   scanf ("%f", &p0);
				   printf ("Digite um novo valor para q0:\n");
				   scanf ("%f", &q0);
			   }
			   
			   k++;
		 }
		 
		 // achar duas raizes (x^2 + p*x + q = 0)
		 // fatorar a eq (m = m-2)
		 if (k<=nmi) {
			printf ("p = %f, q = %f\n", p, q);					// imprimir p e q
			printf ("Resolvendo x^2 + %f x + %f = 0 encontramos x[%d] e x[%d]\n", p, q, j+1, j+2);
			eq2grau (p, q, pr, pi, &j);
			m = m-2;
			for (i=0; i<=m; i++)
				a[m-i] = b[m+2-i];
			printf ("Nova equacao (grau %d):\n", m);			  // imprimir nova eq
			imprimir_eq (a, m);
			if (m>2) {
				if (v[m-2]!=0) {
				p0 = v[m]/v[m-2];
				q0 = v[m-1]/v[m-2];
				}
				else {								  // erro: nao calcula p0 e q0; divisao por zero
					printf ("Nao foi possivel calcular os valores iniciais p0 e q0: divisao por zero\n");
					printf ("Digite um valor para p0:\n");
					scanf ("%f", &p0);
					printf ("Digite um valor para q0:\n");
					scanf ("%f", &q0);
				}
			}
		 }
	}
	
	// se grau = 2 ou 1, entao calcular raizes restantes	
	if (k <= nmi) {
	   if (m==2){
		  p = b[3];
		  q = b[2];
		  eq2grau (p, q, pr, pi, &j);
		  printf ("Resolvendo, encontramos x[%d] e x[%d]\n", n-1, n);
	   }
	   else {
			pr[n-1] = -b[2];
			pi[n-1] = 0;
			printf ("Resolvendo, encontramos x[%d]\n", n);
	   }

	   // imprimir raizes
	   for (i=0; i<=n-1; i++) {
		   printf ("Raiz x[%d] = ", i+1);
		   if (pr[i]!=0 && pi[i]!=0){
			  if (pi[i]>0)
				 printf ("%f + %f i\n", pr[i], pi[i]);
			  else
				 printf ("%f - %f i\n", pr[i], fabs(pi[i]));
		   }
		   else if (pr[i]!=0 && pi[i]==0)
			  printf ("%f\n", pr[i]);
		   else if (pr[i]==0 && pi[i]!=0)
			  printf ("%f i\n", pi[i]);
		   else
			  printf ("%f\n", 0);
	   }
	}
	// erro nmi
	if (k > nmi)
	   printf ("Ultrapassou o numero maximo de iteracoes\n");

}

Arquivado

Equações

Recommended Posts

VictorCacciari 42

Compartilhar este post

Link para o post

Compartilhar em outros sites

Rarylson Freitas 0

Compartilhar este post

Link para o post

Compartilhar em outros sites

VictorCacciari 42

Compartilhar este post

Link para o post

Compartilhar em outros sites

viciado 1

Compartilhar este post

Link para o post

Compartilhar em outros sites

VictorCacciari 42

Compartilhar este post

Link para o post

Compartilhar em outros sites

Rarylson Freitas 0

Compartilhar este post

Link para o post

Compartilhar em outros sites

VictorCacciari 42

Compartilhar este post

Link para o post

Compartilhar em outros sites

Rarylson Freitas 0

Compartilhar este post

Link para o post

Compartilhar em outros sites

Fixar div até atingir uma certa altura

Ajuda com Extends e envio para o banco

PHP+Codeginiter - Orientação para Impressão

Fóruns

Tempo Real

Informação importante