Hello,
I currently have working code to solve a single first order differential equation using a predictor-corrector. I need to modify this to solve a system of first order differential equations.
Here is the system I need to solve:
y'' = -y' + 6y; y(0)=1; y'(0)=-2, on [0,4]
And here is my code to solve a single first order diff. equation (using dev c++)
#include <stdio.h>
#include <vector>
#include <math.h>
#include <cstdlib>
#define MAXDATA 100
using namespace std;
double h, tx[MAXDATA],ty[MAXDATA], x,xi,xf,yi;
int fi,i,n;
double f(double x, double y) {
return (4*x/y - x*y);
}
//classical Runge-Kutta method of order 4
void runge_kutta(double *y) {
double a,b,c,d;
a=h*f(x,*y);
b=h*f(x+h/2,*y+a/2);
c=h*f(x+h/2,*y+b/2);
x+=h;
d=h*f(x,*y+c);
*y =*y + (a+b+b+c+c+d)/6;
}
void equadiff_pc(double *tx,double *ty,double xi,double xf,
double yi,int m,int fi) {
double z=yi,y,w,p[4];
int i,j,k; long ni;
if ((m>MAXDATA) || (fi<1)) return;
h = (xf-xi)/fi/m;
p[3]=f(xi,yi);
tx[0]=xi; ty[0]=yi;
for (k=0, i=1; i<=m; i++) {
ni=(long) (i-1)*fi-1;
for (j=1; j<=fi; j++) {
x=xi+h*(ni+j);
if (++k<4) {
runge_kutta(&z);
p[3-k]=f(x,z);
}
else {
x+=h;
w=z+h/24*(55*p[0]-59*p[1]+37*p[2]-9*p[3]);
do {
y=w;
w=z+h/24*(9*f(x,y)+19*p[0]-5*p[1]+p[2]);
} while (fabs(y-w) > 1e-10);
z=w; p[3]=p[2]; p[2]=p[1];
p[1]=p[0]; p[0]=f(x,z);
}
} // j loop
tx[i]=x; ty[i]=z;
} // i loop
}
int main(void) {
printf("\n Input x begin: "); scanf("%lf",&xi);
printf("\n Input x end : "); scanf("%lf",&xf);
printf("\n Input y begin: "); scanf("%lf",&yi);
printf("\n Input num of points: "); scanf("%d",&n);
printf("\n Input finesse: "); scanf("%d",&fi);
equadiff_pc(tx,ty,xi,xf,yi,n,fi);
//print results
printf("\n");
printf(" X Y \n");
printf(" -----------------------------\n");
for (i=1; i<=n; i++)
printf("%12.6f %12.6f\n",tx[i],ty[i]);
printf(" -----------------------------\n");
system("PAUSE");
return EXIT_SUCCESS;
}