Runge Kutta in C amateur

No momento tento implementar o método de Runge-Kutta de 4º Ordem para equações diferenciais, é um método meio rustico ainda, não está bem lapidado com uma escrita em C, bem organizada e legível. Uma justificativa é o fato de ser iniciante nesta linguagem.
Uma descrição do método de Runge- Kutta de 4º ordem, consiste em determinar constantes apropriadas tais que a fórmula:

yn+1=yn+ak1+bk2+ck3dk4

Que tenha uma certa concordância com um desenvolvimento de Taylor até h⁴ , isto é até o quinto termo. O resultado do método é dado pelas seguintes equações.

yn+1=yn+(k1+k2+k3+k4)/6

k1=h f(xn,yn)

k2=h f(xn+h/2, yn+k1/2)

k3=h f(xn+h/2, yn+k2/h)

k4=h f(xn+h, yn+k3)

Onde h é o passo dado no intervalo de integração, ou seja, se meu intervalo é [0;1] então com um h=0,2, o intervalo terá a seguinte sequência numérica {0, 0,2, 0,4, 0,6, 0,8, 1,0}.
O código fonte em C que trago aqui é para uma equação em específico. Ainda estou trabalhando em um código que eu somente coloque a função sem ter que modificar o código fonte.

Segue o código em C do Método de Runge Kutta, para a resolução da equação:

dy/dx = x-y (1)
ou
y'=x-y

/*********************************************************************/
/*Bibliotecas*/
#include
#include
#include
FILE*rk4;
/*Inicio do programa*/
main(){
rk4=fopen("rk4.txt","w");
int n;
double x,yv,h,yn,x0;
double f1,f2,f3,f4;

h=0.2;
x0=0;
yv=2;

for(n=0;n<=10;n++){
/*Variaveis do método dy/dx=x+y, aqui já na equação diferencial a ser integrada*/

f1=h*(x+yv);
f2=h*((x+h/2)+(yv+f1/2));
f3=h*((x+h/2)+(yv+f2/2));
f4=h*((x+h)+(yv+f3));

/* Escreve no terminal e no arquivo*/

fprintf(rk4,"%lf%18lf\n",x,yn);
printf("%lf%18lf\n",x,yn); }

/**** Escreve no terminal e no arquivo*****/

fprintf(rk4,"%lf%18lf\n",x,yn); printf("%lf%18lf\n",x,yn);

}
fclose(rk4); }

/*Fim do programa*/
/************************************************************************/

A solução algébrica da equação (1)

y(x)=3exp(x)-x-1 (2)

A seguir os gráficos da equação resolvida pelo método e algebricamente, apesar dos erros obtidos pelos calculos computacionais, podemos notar que o seu bom resultado gráficamente.

Gráfico calculado pelo método de Runge Kutta de quarta ordem














Gráfico obtido pela função exata calculada algebricamente.

Programa (2)

O código fonte aqui descrito neste post, é uma rotina para cálculos de volumes de esferas. Um código básico mas que pode ser util para quem esta aprendendo a implementar cálculos em linguagem C.

/*********Bibliotecas********/
#include
#include
/**************************/
/************Inicio do programa**************/
main(){
double V,r,pi,Vf,d;

pi=3.14159;
printf("Digite o raio da esfera:\n");
scanf("%lf",&r);
d=2*r;
V=(4/3)*pi*pow(r,3);
printf("Volume da esfera:%lf u.v\n",V);
Vf=pow(d,3)-(pi/6)*pow(d,3);
printf("Agora se imaginar extrair uma esfera de dentro de um cubo,\n o volume da figura:%lf\n",Vf);
return 0;
}
/************Fim do programa********************/

Linguagem C

Depois de tanto pesquisar sobre linguagens de programação optei pela qual o meu grupo de pesquisa trabalha, isso é devido a forças maiores. Explicações de quem a criou, desenvolveu é encontrado em qualquer site de busca como o Google. Aqui pra quem interessar é mais um blog de discussões da linguagem C aplicada a matemática e a física mais precisamente na dinâmica não-linear de sistemas complexos.

Aqui é uma rotina introdutória, com forma o tempo evoluir eu vou colocando mais.
A rotina a seguir é uma rotina que calcula o montante, um exemplo simples para matemática financeira, o compilador utilizado foi o GCC e o sistema operacional é o UBUNTU 10.04 (Lucid Lynx) que pode ser obtidos (aqui).

[rotina]

#include
#include

main(void)
{
int n;
float a,p,r,m;
printf("Digite o capital a ser investido:\n");
scanf("%f",&p);
printf("Digite a taxa que será aplicada:\n");
scanf("%f",&r);
{
printf("%4s%21s%20s\n","Tempo","Juro composto","Montante");
for(n=1;n<=10;n++)
{a=p*pow(1+(r/100),n);
m=p+a;
fprintf("%4d%21f%21f\n",n,a,m);
}
}
return 0;
}