C/C++ Runge Kutta Script | C++ | Beginner

[okkvlt]

Hi. I am working on a runge kutta/integral finder. my goal is to make it "exact to the dx", so to speak.
im a beginner at c++, so bear with me.
the compiler I am using is dev cpp, incase it matters.

my main problem was that apparently i can't have a long double array of size greater than 160,000. this would limit me to an x array of 80k and a y array of 80k. (meaning i could only do 80k iterations of runge kutta). i worked around this by making an x and y array of 40k entries each, and then subdividing each increment into a maximum of 40k entries, giving me the maximum of up to 1.6billion runge kutta iterations. (although 40million works fine, and i don't feel like waiting several hours for my POS 1.3ghz computer to do 1.6 billion runge kutta iterations). the only drawback is that after each completion of an interval the fine-precision of the subintervals is lost, i worked around this for the print f[x] and definite integral part, but its kind of messy.

Some questions:
how do i use scanf to declare the function to be integrated/solved? it would be real nice if i could do this in the command prompt window, instead of having to go into the compiler and edit the program.

how do i fprintf or printf a long double? it seems as though the print function only works for double precision, not long double, so i have to convert my long double to merely double precision before i can see it or save it :( .

How do i scanf a double or long double? scanf only seems to work for float numbers; this seriously affects things i want to do with the program, like suppose i want to find the definite integral over some multiple of pi. then id have to go into the compiler and edit the program.

how do i use the second derivative to optimize numerical analysis? i know that if the second derivative is large then i should use a smaller increment size, but i can't figure out how to implement this into the program.

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


/*define f[x]*/
long double f(long double x, long double y)
{
long double f=cos(atan(y/x))*pow(pow(x,2)+pow(y,2),.5);
return f;
}
int main()
{
/*get values*/
long int num;
float xinit, yinit, xend, xa,xb;
long double x[40000], y[40000], h;
printf("find y[x] given dy/dx=f[x,y]\n");
printf("initial x=");
scanf("%f", &xinit);
printf("initial y=");
scanf("%f", &yinit);
printf("final x=");
scanf("%f", &xend);
printf("number of subincrements=");
scanf("%d", &num);
double dxdisp=(xend-xinit)/(40000*num), xdisp, ydisp;
printf("dx=%.13f\n", dxdisp);
x[0]=xinit;
y[0]=yinit;
long int i=0, n=0;
long double xsub[num],ysub[num],a,b,c,d, dx=(xend-xinit)/(40000*num), yb, ya;
while(i<=40000)
{
n=0;
ysub[0]=y;
xsub[0]=x;
while(n<=num)
{
a=dx*f(xsub[n],ysub[n]);
b=dx*f(xsub[n]+.5*dx,ysub[n]+.5*a);
c=dx*f(xsub[n]+.5*dx,ysub[n]+.5*b);
d=dx*f(xsub[n]+dx,ysub[n]+c);
ysub[n+1]=ysub[n]+(a+2*b+2*c+d)/6;
xsub[n+1]=xsub[n]+dx;
n=n+1;
}
x[i+1]=xsub[num];
y[i+1]=ysub[num];
i=i+1;
printf("%d/40000 completed\r", i);
}
printf("-------------------------------\n");
long int choose=10;
long int display;

printf(" 0:exit\n 1:save data\n 2:print y[x]\n 3:definite integral if f=f[x]\n");

while(choose!=0)
{
choose=10;
printf("choice=");
scanf("%d", &choose);

/*save*/
if(choose==1)
{
display=0;
FILE *fp;
fp=fopen("graph.txt","w");
while(display<=40000)
{
xdisp=x[display];
ydisp=y[display];
fprintf(fp, "%.13f %.13f\n", xdisp, ydisp);
display=display+4;
}
int fclose(FILE *fp);
printf("finished saving\n");
}
/*save*/

/*y[x]*/
if(choose==2)
{
printf("x=");
scanf("%f", &xa);
i=floor((xa-xinit)/(dx*num));
h=(xa-x)/(num);
n=0;
xsub[0]=x;
ysub[0]=y;
while(n<=num)
{
a=h*f(xsub[n],ysub[n]);
b=h*f(xsub[n]+.5*h,ysub[n]+.5*a);
c=h*f(xsub[n]+.5*h,ysub[n]+.5*b);
d=h*f(xsub[n]+h,ysub[n]+c);
ysub[n+1]=ysub[n]+(a+2*b+2*c+d)/6;
xsub[n+1]=xsub[n]+h;
n=n+1;
}
xdisp=xsub[num];
ydisp=ysub[num];
printf("y(%.13f)=%.13f\n\n", xdisp, ydisp);
}
/*defint*/
if(choose==3)
{
printf("b=");
scanf("%f", &xb);
printf("a=");
scanf("%f", &xa);
i=floor((xa-xinit)/(dx*num));
h=(xa-x)/(num);
n=0;
xsub[0]=x;
ysub[0]=y;
while(n<=num)
{
a=h*f(xsub[n],ysub[n]);
b=h*f(xsub[n]+.5*h,ysub[n]+.5*a);
c=h*f(xsub[n]+.5*h,ysub[n]+.5*b);
d=h*f(xsub[n]+h,ysub[n]+c);
ysub[n+1]=ysub[n]+(a+2*b+2*c+d)/6;
xsub[n+1]=xsub[n]+h;
n=n+1;
}
ya=ysub[num];
i=floor((xb-xinit)/(dx*num));
h=(xb-x)/(num);
n=0;
xsub[0]=x;
ysub[0]=y;
while(n<=num)
{
a=h*f(xsub[n],ysub[n]);
b=h*f(xsub[n]+.5*h,ysub[n]+.5*a);
c=h*f(xsub[n]+.5*h,ysub[n]+.5*b);
d=h*f(xsub[n]+h,ysub[n]+c);
ysub[n+1]=ysub[n]+(a+2*b+2*c+d)/6;
xsub[n+1]=xsub[n]+h;
n=n+1;
}
yb=ysub[num];
ydisp=yb-ya;
printf("%f\nSf[x]dx=%.13f\n%f\n\n", xb,ydisp,xa);
}
/*defint*/
}
/*menu*/


/*end*/
}

 

[CFDFEAGURU]

First you are attempting the impossible, you can only be "exact" to the simplest of systems. Real world problems can only be approximated (sometimes to a very high degree) and even then you can only be as accurate as the machine precision in your CPU. Secondly, RK4 is never going to get you high precision with a real equation or system of equations. You need to use a much more advanced algorithm like 8th order Runge-Kutta. Also, to achieve any hope of fast convergence you need an automated time stepper in your algorithm.

Lastly, you are doing something that has already been done for you. Go to the Numerical Recipe site. www.nr.com and read their methods of numerical solutions to ODE's. Also, you aren't losing programming knowledge by using their routines. In fact, you will gain a lot by just following them and seeing how they work.

Thanks
Matt

 

[rcgldr]

RK4 is never going to get you high precision with a real equation or system of equations. You need to use a much more advanced algorithm like 8th order Runge-Kutta.

An alternative to higher order algorithms is to use a smaller step size. Also, consider using an implicit (predictor-corrector) algorithm, which uses the next step value as an input, by making an initial guess for the next step value, then iterating using sucessive "guesses" from previous steps.

I explained a simple implicit algorithm here:

https://www.physicsforums.com/showpost.php?p=2285015&postcount=10

Another issue is that summation of a large number of values results in truncation, although it woud seem that long doubles (80 bit?) numbers would take care of this. There are algorithms that track truncation errors, similar to: truncation = (number-((sum+number)-sum))). Another option is to use a summation function. One implementation is to use an array of floating point numbers indexed by the exponent of the number. The array is initialized to zeros, and the function called for each new number to be added to the array. The exponent of the new number is used to index into the array. If the entry is zero, the new number is stored. If the entry is non-zero, the new number is added to the array entry, the array entry is zeroed, and the process is repeated using the sum's exponent to index the array until a zero entry is found or overflow occurs. There is a another function that returns the current sum by adding the numbers in array (exponent) order.