Código em C não dando o resultado esperado

[LucianoPeixoto_rj 0]

Bom dia,

 

Sou novo em c e não estou conseguindo entender o que há de errado no meu código de Runge-Kutta 4, fiz ele em c e em pascal, o código é o mesmo só adaptei a sintaxe, entretanto os resultados dão diferentes. Os resultados do código em pascal estão certos porém o do C não. Dica pra quem for me ajudar, só verifique as três colunas do último resultado, se derem iguais está tudo certo. Obrigado!

 

>>>>>>>CÓDIGO EM C<<<<<<<

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double f(double t, double x, double v, int op)
{   switch (op)
    {
       case 1 : 
            return (0*t)+(0*x)+(1*v); 
            break; 
       case 2 : 
            return ((24/20)*sin(15.0*t))-((20*20)*x)-(2*20*0.1625*v); 
            break;
    }
}
#define nper 2001
#define nvar 3
#define neq 2
#define incr 0.001
int main(int argc,  char** argv) 
{   FILE *fp = fopen("file.txt", "w"); 
    long int n, m, eq;
    double x[nvar], k1[neq], k2[neq], k3[neq], k4[neq];
    double y[nper+1][nvar];
    y[1][1] = 0.0; 
    y[1][2] = 0.0; 
    y[1][3] = 0.1; 
    for (n=1; n<=nper; ++n)
    {    
        for (eq=1; eq<=neq; ++eq)
        {
            for (m=1; m<=nvar; ++m)            
                x[m]=y[n][m];           
            k1[eq] = f(x[1],x[2],x[3],eq);
        }
        eq=1; m =1;
        for (eq=1; eq<=neq; ++eq)
        {
            for (m=1; m<=nvar; ++m)
            {  
                if (m==1)       
                   x[m] = y[n][m] + 0.5*incr;
                else
                   x[m] = y[n][m] + 0.5*incr*k1[m-1];
            }   
            k2[eq] = f(x[1],x[2],x[3],eq);
        }
        eq=1; m =1;
        for (eq=1; eq<=neq; ++eq)
        {
            for (m=1; m<=nvar; ++m)
            {
                if (m==1)       
                   x[m] = y[n][m] + 0.5*incr;
                else
                   x[m] = y[n][m] + 0.5*incr*k2[m-1];                    
            }
            k3[eq] = f(x[1],x[2],x[3],eq);
        }
        eq=1; m =1;
        for (eq=1; eq<=neq; ++eq)
        {
            for (m=1; m<=nvar; ++m)
            {
                if (m==1)       
                   x[m] = y[n][m] + incr;
                else
                   x[m] = y[n][m] + incr*k3[m-1];
            }
            k4[eq] = f(x[1],x[2],x[3],eq);        
        } 
        m=1;  
        for (m=1; m<=nvar; ++m)
        {
         if (m==1)       
            y[n+1][m] = y[n][m] + incr;
         else
            y[n+1][m] = y[n][m]+(incr/6)*(k1[m-1]+2*k2[m-1]+2*k3[m-1]+k4[m-1]);
        } 
    }
    for (n=1; n<=nper; ++n)
        {printf("%g, %g, %g\n", y[n][1], y[n][2], y[n][3]);
         fprintf(fp, "%g    %g    %g\n", y[n][1], y[n][2], y[n][3]);
        }
    fclose(fp);
    system("pause");
    return (EXIT_SUCCESS);        
} 

>>>>>>>CÓDIGO EM PASCAL<<<<<<<

Program Runge_Kutta;
Const
	nper = 2001;
	nvar = 3;
	neq = 2;
	incr = 0.001;
Var
	y: array [1..2002,1..nvar] of real;
	x: array [1..nvar] of real;
	k1,k2,k3,k4: array[1..neq] of real;
	n,m,eq:integer;
	s: Text;
Function f(t,x,v:real;c:integer):real; 
	Begin
		if c = 1 then
			f:=(0*t)+(0*x)+(1*v)				
		else
   			f:=((24/20)*sin(15*t))-((20*20)*x)-(2*20*0.1625*v);
	End;
 Begin
 	assign(s, 'resultado1.txt');
 	rewrite(s);
 	y[1,1]:=0;
 	y[1,2]:=0;
 	y[1,3]:=0.1;
 	for n:=1 to nper do
 	begin
 		for eq:=1 to neq do
 		begin
 				for m:=1 to nvar do
 					x[m]:= y[n,m];
 				k1[eq]:= f(x[1],x[2],x[3],eq)
 		end;
 		for eq:=1 to neq do
 		begin
 				for m:=1 to nvar do
 				begin
 					if m=1 then
 						x[m]:= y[n,m] + 0.5*incr
					else
						x[m]:= y[n,m] + 0.5*incr*k1[m-1];	
 				end;	
 				k2[eq]:= f(x[1],x[2],x[3],eq);
 		end;
 		for eq:=1 to neq do
 		begin
 				for m:=1 to nvar do
 				begin
 					if m=1 then
 						x[m]:= y[n,m] + 0.5*incr
					else
						x[m]:= y[n,m] + 0.5*incr*k2[m-1];	
 				end;	
 				k3[eq]:= f(x[1],x[2],x[3],eq);
 		end; 		
 		for eq:=1 to neq do
 		begin
 				for m:=1 to nvar do
 				begin
 					if m=1 then
 						x[m]:= y[n,m] + incr
					else
						x[m]:= y[n,m] + incr*k3[m-1];	
 				end;	
 				k4[eq]:= f(x[1],x[2],x[3],eq);
 		end; 		
		for m:=1 to nvar do
		begin
			if m = 1 then
				y[n+1,m]:= y[n,m]+ incr
			else
				y[n+1,m]:= y[n,m]+(incr/6)*(k1[m-1]+2*k2[m-1]+2*k3[m-1]+k4[m-1]);
		end;
	end;
		for n:=1 to nper do
		begin
			writeln(y[n,1],' ',y[n,2],' ',y[n,3]);
			writeln(s,y[n,1],' ',y[n,2],' ',y[n,3]);
		end;

 End.

 

[_Isis_ 202]

No Pascal, por acaso, as variáveis não inicializadas assumem valores default?

[LucianoPeixoto_rj 0]

Sim, assumem valor 0.

[LazaroBinda 3]

Qual ide vc está usando pra pascal??? Tenho minhas dúvidas se ele carrega valor default.

[LucianoPeixoto_rj 0]

Lazarus Free Pascal, mas o que isso implica no código em C ?

[_Isis_ 202]

Porque em C não existe inicialização automática de variáveis. Quando vc declara uma variável tem um valor aleatório na posição de memória. Se o teu algoritmo precisa dos valores iniciais iguais a zero, vc vai ter que atribuir manualmente.

[LucianoPeixoto_rj 0]

Entendi, eu tinha percebido por isso eu coloquei entre dois for (eq=1; m =1;) mas não teve efeito nenhum o resultado foi o mesmo. Do código as variáveis que precisavam ter um valor antes de entrar no for eu declarei, que foram y[1][1], y[1][2] e y[1][3], já as outras são inicializadas dentro do for quando coloco for eq = 1 ... por exemplo, então não vejo nenhuma variável não declarada, mas obrigado e continuem tentando porque ta ruim de sair isso. rsrsrs

[_Isis_ 202]

Variável não declarada não é a mesma coisa que variável não inicializada.

Provavelmente é diferença na precisão do ponto flutuante, porque a coisa já começa a ficar diferente aqui:

for eq:=1 to neq do
        begin
                for m:=1 to nvar do
                    x[m]:= y[n,m];
                k1[eq]:= f(x[1],x[2],x[3],eq)
        end;

 

Na saída do pascal o formato é 0.00000000000E+000, enquanto que na saída do C vc enxerga um 4.10786e-39 (o que seria basicamente igual a zero).

Problema é que não consigo nem encaixar minimamente teu código no método por causa dos nomes das variáveis e porque não há enunciado... http://pt.wikipedia.org/wiki/M%C3%A9todo_de_Runge-Kutta

[LucianoPeixoto_rj 0]

Hum.. é diferente mas não acho que isso seja a causa por ser um valor extremamente próximo de zero. Esse método que você postou é válido para uma equação apenas, e o meu caso é um sistema. É parecido, mas tem um pulo do gato rs. Eu vou escrever o método para o sistema com as variáveis que eu usei no código e vou explicar aonde as linhas do código entram no método. Isso vai dar uma baita ajuda, já que está difícil de achar o erro, eu pensava que tinha errado algo estúpido que na primeira resposta já ia aparecer o erro, mas não está sendo o caso, então mais tarde postarei o método.

[LucianoPeixoto_rj 0]

Agradeço o interesse e a ajuda de todos, já resolvi o problema, primeiro botei os vetores para começar em 0 pra ficar mais apropriado ao c e ajustei os limites deles, ai percebi que não estava pegando os "Ks" certos e acho q isso estava acontecendo no anterior também, então ajustei a posição dos "Ks" para [m-1] e pronto, os resultados do C e do Pascal agora batem.

 

Com relação ao que falei sobre o que vc postou Isis aquele é para uma equação só, caso precise resolver um sistema de equações, cada uma delas terá seu conjunto de K1,K2,K3 e K4 porém devem ser intercambiáveis conforme o link

 

https://www.nsc.liu.se/~boein/f77to90/rk.html

CÓDIGO CORRETO

/* 
 * File:   main.c
 * Author: Luciano
 *
 * Created on 13 de Outubro de 2014, 11:44
 */

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
double f(double t, double x, double v, int op)
{   switch (op)
    {
       case 0 : 
            return (0*t)+(0*x)+(1*v); 
            break; 
       case 1 : 
            return (1.2*sin(15*t))-((20*20)*x)-(2*20*0.1625*v); 
            break;
    }
}
#define nper 2000
#define nvar 3
#define neq 2
#define incr 0.001
int main(int argc,  char** argv) 
{   FILE *fp = fopen("file.txt", "w"); 
    long int n, m, eq;
    double x[nvar], k1[neq], k2[neq], k3[neq], k4[neq];
    double y[nper+1][nvar];
    y[0][0] = 0.0; 
    y[0][1] = 0.0; 
    y[0][2] = 0.1; 
    n = 0;
    eq = 0;
    m = 0;
    for (n=0; n<=nper; n=n+1)
    {
        printf("%g, %g, %g\n", y[n][0], y[n][1], y[n][2]);
        fprintf(fp, "%g    %g   %g\n", y[n][0], y[n][1], y[n][2]);    
        for (eq=0; eq<=(neq-1); eq=eq+1)
        {
            for (m=0; m<=(nvar-1); m=m+1)       
                x[m]=y[n][m];   
            k1[eq] = f(x[0],x[1],x[2],eq);
        }
        for (eq=0; eq<=(neq-1); eq=eq+1)
        {
            for (m=0; m<=(nvar-1); m=m+1)
            {  
                if (m==0)       
                   x[m] = y[n][m] + 0.5*incr;
                else
                   x[m] = y[n][m] + 0.5*incr*k1[m-1];
            }   
            k2[eq] = f(x[0],x[1],x[2],eq);
        }
        for (eq=0; eq<=(neq-1); eq=eq+1)
        {
            for (m=0; m<=(nvar-1); m=m+1)
            {
                if (m==0)       
                   x[m] = y[n][m] + 0.5*incr;
                else
                   x[m] = y[n][m] + 0.5*incr*k2[m-1];                    
            }
            k3[eq] = f(x[0],x[1],x[2],eq);
        }
        for (eq=0; eq<=(neq-1); eq=eq+1)
        {
            for (m=0; m<=(nvar-1); m=m+1)
            {
                if (m==0)       
                   x[m] = y[n][m] + incr;
                else
                   x[m] = y[n][m] + incr*k3[m-1];
            }
            k4[eq] = f(x[0],x[1],x[2],eq);        
        }   
        for (m=0; m<=(nvar-1); m=m+1)
        {
            if (m==0)       
               y[n+1][m] = y[n][m] + incr;
            else
            y[n+1][m] = y[n][m]+(incr/6)*(k1[m-1]+2*k2[m-1]+2*k3[m-1]+k4[m-1]);
        } 
    }
    fclose(fp);
    system("pause");
    return (EXIT_SUCCESS);   
     
}