- #1
okkvlt
- 53
- 0
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*/
}
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*/
}