Hi. im working on a runge kutta/integral finder. my goal is to make it "exact to the dx", so to speak.(adsbygoogle = window.adsbygoogle || []).push({});

im a beginner at c++, so bear with me.

the compiler im using is dev cpp, incase it matters.

my main problem was that apparently i cant 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 dont 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 cant 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*/

}

**Physics Forums - The Fusion of Science and Community**

The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

# Runge kutta script

Loading...

Similar Threads - Runge kutta script | Date |
---|---|

Fortran Runge - Kutta 4 understanding | Oct 17, 2017 |

Python Runge Kutta for nonlinear system of equation | Oct 3, 2016 |

Plotting Orbits on MATLAB | Aug 16, 2016 |

C/++/# Runge-Kutta 4th order program copying array element | Aug 9, 2016 |

C/++/# Problems with Runge-Kutta method | Jul 24, 2016 |

**Physics Forums - The Fusion of Science and Community**