PDA

View Full Version : Runge-Kutta in C++


Pauly Man
Apr23-03, 02:39 AM
Hi guys.

I was wondering if anyone could help me with this code. I have sucessfully created a program in visual basic that can run a runge-kutta method. However I want to create one in c++, maybe eventually turn it into a .dll when i work out how to create and use them. [:)]

Here is the code:


// Runge-Kutta.cpp
//--------------------------------------------------
//A Runge-Kutta Method for solving Differential Equations
//of the form y'=f(x,y) ; y(x0)=y0
//--------------------------------------------------

#include <iostream>
#include <iomanip>
using namespace std;

//Define constants
#define X0 0
#define Y0 0
#define H 0.2
#define N 5

//Define Functions
double f(double x, double y);
double runge(double x, double y);


//Main Function
int main(double x, double y)
{
cout<<"\t*** Euler Method ***"
<<"\n\n";
cout<<" "
<<setw(12)<<"x"<<setw(12)<<"\ty"
<<"\n"
<<"\t------------------------------"
<<"\n";
y=Y0;
for(int i=0;i<=5;i++)
{
x=X0+(i*H);
y=runge(x,y);
cout<<left<<setw(6)<<i<<"|"
<<setprecision(4)<<left<<setw(8)<<"\t"<<x
<<setprecision(4)<<left<<setw(8)<<"\t"<<y;
cout<<"\n\n";
}
cout<<"\n\n";
return 0;
}

double runge(double x, double y)
{
double K1 = (H * f(x,y));
double K2 = (H * f((x + 1 / 2 * H), (y + 1 / 2 * K1)));
double K3 = (H * f((x + 1 / 2 * H), (y + 1 / 2 * K2)));
double K4 = (H * f((x + H), (y + K3)));
double runge = (y + (1 / 6) * (K1 + 2 * K2 + 2 * K3 + K4));
return runge;
}

double f(double x, double y)
{
double f = x+y;
return f;
}


It is supposed to print out a table, with the x values in one column, and y values in the other. Each successive iteration is in a new row of the table.

So far so good, the code creates the table fine.

I created a similiar program to run a Euler method, and it works just fine. But in the Runge-Kutta it just prints the y values as zero. I must have done something wrong somewhere.

Pauly Man
Apr23-03, 03:56 AM
I've worked it out now.

It was the final line of the "runge" function.

I had the brackets wrong, it should read:

double runge = (y + ((K1 + 2 * K2 + 2 * K3 + K4)/6));

jonnylane
Apr23-03, 09:09 AM
i wrote a very similar program in FORTRAN 90 a couple of years ago.

just out of interest, here's the subroutine for the program it was used in:

SUBROUTINE RUNKUT(H,F,X,N,Y)
! Performs one step of a fourth order Runge-Kutta integration
IMPLICIT NONE
INTEGER, PARAMETER :: DOUB_PREC=SELECTED_REAL_KIND(P=10,R=30)
INTEGER :: N ! number of y's
REAL(DOUB_PREC) :: H,X,Y(N) ! x step size, x, y's
EXTERNAL F ! subroutine of the form F(X,Y,YPRIME)
! The subroutine F must return derivatives of Y's wrt X in YPRIME
! Both the independent variable X and the dependent variables Y are
! updated on each call.
! Local automatic arrays for workspace
REAL(DOUB_PREC) :: F1(N),F2(N),F3(N),F4(N)
REAL(DOUB_PREC) :: Y1(N),Y2(N),Y3(N),Y4(N)
REAL(DOUB_PREC) :: K1(N),K2(N),K3(N),K4(N)
REAL(DOUB_PREC) :: X1,X2,X3,X4
! Define steps for 4th order terms
X1=X
X2=X1+H*0.5
X3=X2
X4=X1+H
! Calculate 4th order terms
Y1=Y
CALL F(X1,Y1,F1)
K1=H*F1
Y2=Y1+K1*0.5
CALL F(X2,Y2,F2)
K2=H*F2
Y3=Y1+K2*0.5
CALL F(X3,Y3,F3)
K3=H*F3
Y4=Y1+K3
CALL F(X4,Y4,F4)
K4=H*F4
! Perform 4th order step
Y=Y1+(K1+2.0D0*K2+2.0D0*K3+K4)/6.0D0
! Update independent variable
X=X4
END SUBROUTINE RUNKUT

acemv
May1-03, 09:50 PM
you're all set, i think you should overload some functions though.

efebest
Nov8-08, 02:27 PM
hey dude nice codes there but can i get ur codes in vb am working on something similar thnx

D H
Nov8-08, 02:49 PM
This thread is over five years old! :surprised

None of the people who posted in this thread are around to give you their code.