C/C++ How Can I Fix My C++ Runge-Kutta Method to Correctly Update Y Values?

[Pauly Man]

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 similar 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]

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]

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]

you're all set, i think you should overload some functions though.

 

[efebest]

hey dude nice codes there but can i get ur codes in vb am working on something similar thnx

 

[D H]

This thread is over five years old!

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