Need help with a research project

Page 1 of 1

0 Replies - 1000 Views - Last Post: 10 December 2016 - 09:03 PM

#1 hiephoang  Icon User is offline

  • New D.I.C Head

Reputation: -1
  • View blog
  • Posts: 23
  • Joined: 18-December 15

Need help with a research project

Posted 10 December 2016 - 09:03 PM

I'm doing research concerning Heat and Mass Transfer.

The program created in order to solve a system of 3 Ordinary Differential Equations using Runge Kutta and shooting method:

S'''+3*S*S''-2*(S'*S')+P-0.5*F=0
P''+3*10*S*P'+0.5*F'*P'+0.5(P'*P')=0
F''+3*10*10*S*F'+P''=0

Initial and boundary conditions:
x=0: S=0, S'=0, P=1, F=1
x= infinity: S'=0, P=0, F=0
Is there any one who can help me with this program? The results produced contain some errors as following:

number of interval in the coordinate n= 300000
step in the coordinate h= 3.33333e-005
bb= 2.84057
cc= 0.777512
ee= -3.45846e-323
iter = 1 dd= 3.61808
bb= -1.#IND
cc= -1.#IND
ee= -1.#IND
iter = 2 dd= 1.#QNAN
x[i] = 0then u[i]= 0
x[i] = 0then y[i]= 0
x[i] = 0then z[i]= -1.#IND
x[i] = 0then w[i]= 1
x[i] = 0then p[i]= -1.#IND
x[i] = 0then m[i]= 1
x[i] = 0then v[i]= -1.#IND
x[i] = 3.33333e-005then u[i]= -1.#IND
x[i] = 3.33333e-005then y[i]= -1.#IND
x[i] = 3.33333e-005then z[i]= -1.#IND
x[i] = 3.33333e-005then w[i]= -1.#IND
x[i] = 3.33333e-005then p[i]= -1.#IND
x[i] = 3.33333e-005then m[i]= -1.#IND
x[i] = 3.33333e-005then v[i]= -1.#IND
x[i] = 6.66667e-005then u[i]= -1.#IND
x[i] = 6.66667e-005then y[i]= -1.#IND
x[i] = 6.66667e-005then z[i]= -1.#IND
x[i] = 6.66667e-005then w[i]= -1.#IND
x[i] = 6.66667e-005then p[i]= -1.#IND
x[i] = 6.66667e-005then m[i]= -1.#IND
x[i] = 6.66667e-005then v[i]= -1.#IND
x[i] = 0.0001then u[i]= -1.#IND
x[i] = 0.0001then y[i]= -1.#IND
x[i] = 0.0001then z[i]= -1.#IND
x[i] = 0.0001then w[i]= -1.#IND
x[i] = 0.0001then p[i]= -1.#IND
x[i] = 0.0001then m[i]= -1.#IND
x[i] = 0.0001then v[i]= -1.#IND
x[i] = 0.000133333then u[i]= -1.#IND
x[i] = 0.000133333then y[i]= -1.#IND
x[i] = 0.000133333then z[i]= -1.#IND
x[i] = 0.000133333then w[i]= -1.#IND
x[i] = 0.000133333then p[i]= -1.#IND
x[i] = 0.000133333then m[i]= -1.#IND
x[i] = 0.000133333then v[i]= -1.#IND 


The code is following:
#include<math.h>
double f1(double e1, double e2, double e3, double e4, double e5, double e6, double e7, double e8)
{
	return e3;
}
//---------------------
double f2(double e1, double e2, double e3, double e4, double e5, double e6, double e7, double e8)
{
	return e4;
}
//---------------------
double f3(double e1, double e2, double e3, double e4, double e5, double e6, double e7, double e8)
{
	return  -3.0*e2*e4 + 2.0*e3*e3 - e5 + 0.5*e7;
}
//---------------------
double f4(double e1, double e2, double e3, double e4, double e5, double e6, double e7, double e8)
{
	return e6;
}
double f5(double e1, double e2, double e3, double e4, double e5, double e6, double e7, double e8)
{
	return -30.0*e2*e6-0.5*e8*e6-0.5*e6*e6;
}
double f6(double e1, double e2, double e3, double e4, double e5, double e6, double e7, double e8)
{
	return e8;
}
double f7(double e1, double e2, double e3, double e4, double e5, double e6, double e7, double e8)
{
	return -3.0*10*10*e2*e8+3.0*10*e2*e6+0.5*e8*e6+0.5*e6*e6;
}
#include<stdio.h>
#include<conio.h>
#include<iostream>
#include<fstream>
using namespace std;

//------------------------
int main(void)
{
	ofstream scientificwork2("scientificwork2.txt");
	const double EPS =1e-10;
	int i, n = 300000;
	double h=10.0/n;
	scientificwork2 << "number of interval in the coordinate n= " << n << endl;
	//
	//
	//
	double *x;
	x = new double[n + 1];
	x[0] = 0.0;
	//
	for (i = 0; i<n + 1; i++)
		x[i] = i*h;
	//
	//
	scientificwork2<<"step in the coordinate h= "<<h<<endl;
	//
	double *u;
	u = new double[n + 1];
	u[0] = 0.0;
	double *y;
	y = new double[n + 1];
	y[0] = 0.0;
	double *z;
	z = new double[n + 1];
	z[0] = 0.7;
	double *w;
	w = new double[n + 1];
	w[0] = 1.0;
	double *p;
	p = new double[n + 1];
	p[0] = -0.3;
	double *m;
	m = new double[n + 1];
	m[0] = 1.0;
	double *v;
	v = new double[n + 1];
	v[0] = -1.25;
    //
	int iter;
	//
	double eta = 0.7, ep = -0.3, ha = -1.25;
	double hz, hp, hv, dd, dz, dp, dv;
	double bb, cc, ee, bbz, ccz, eez, bbp, ccp, eep, bbv, ccv, eev ;
	double a1, a2, a3, a4;
	double b1, b2, b3, b4;
	double c1, c2, c3, c4;
	double d1, d2, d3, d4;
	double g1, g2, g3, g4;
	double r1, r2, r3, r4;
	double t1, t2, t3, t4;
	iter = 0;
	do
	{
		iter++;
		u[0] = 0.0;
		y[0] = 0.0;
		z[0] = eta;
		w[0] = 1.0;
		p[0] = ep;
		m[0] = 1.0;
		v[0] = ha;
		//
		//
		for (i = 0; i<n; i++)

		{
			a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			r1 = f6(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			t1 = f7(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
            //
			a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			r2 = f6(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			t2 = f7(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			//
			//		f1 returned	0.25000000000000000	double
			//
			a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			r3 = f6(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			t3 = f7(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			//
			//
			a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			r4 = f6(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			t4 = f7(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;

			//
			u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
			y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
			z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
			w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
			p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;
			m[i + 1] = m[i] + (r1 + 2.0*r2 + 2.0*r3 + r4) / 6.0;
			v[i + 1] = v[i] + (t1 + 2.0*t2 + 2.0*t3 + t4) / 6.0;
		}
		//
		bb = y[n] - 0.0;
		cc = w[n] - 0.0;
		ee = v[n] - 0.0;
		scientificwork2 << "bb= " << bb << endl;
		scientificwork2 << "cc= " << cc << endl;
		scientificwork2 << "ee= " << ee << endl;
		//
		//
		u[0] = 0.0;
		y[0] = 0.0;
		if (fabs(eta) <= 1.0) hz = sqrt(EPS);
		if (fabs(eta)>1.0) hz = eta*sqrt(EPS);
		z[0] = eta + hz;
		w[0] = 1.0;
		p[0] = ep;
		m[0] = 1.0;
		v[0] = ha;
		//
		//
			//
		for (i = 0; i<n; i++)

		{
			a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			r1 = f6(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			t1 = f7(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
            //
			a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			r2 = f6(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			t2 = f7(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			//
			//		f1 returned	0.25000000000000000	double
			//
			a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			r3 = f6(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			t3 = f7(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			//
			//
			a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			r4 = f6(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			t4 = f7(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			//
			u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
			y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
			z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
			w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
			p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;
			m[i + 1] = m[i] + (r1 + 2.0*r2 + 2.0*r3 + r4) / 6.0;
			v[i + 1] = v[i] + (t1 + 2.0*t2 + 2.0*t3 + t4) / 6.0;
		}
		bbz = y[n] - 0.0;
		ccz = w[n] - 0.0;
		eez = v[n] - 0.0;
		//
		//
		u[0] = 0.0;
		y[0] = 0.0;
		z[0] = eta;
		w[0] = 1.0;
		if (fabs(ep) <= 1.0) hp = sqrt(EPS);
		if (fabs(ep)>1.0) hp = ep*sqrt(EPS);
		p[0] = ep + hp;
		m[0] = 1.0;
		v[0] = ha;
        //
		//
			for (i = 0; i<n; i++)

		{
			a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			r1 = f6(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			t1 = f7(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
            //
			a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			r2 = f6(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			t2 = f7(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			//
			//		f1 returned	0.25000000000000000	double
			//
			a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			r3 = f6(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			t3 = f7(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			//
			//
			a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			r4 = f6(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			t4 = f7(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			//
			u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
			y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
			z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
			w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
			p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;
			m[i + 1] = m[i] + (r1 + 2.0*r2 + 2.0*r3 + r4) / 6.0;
			v[i + 1] = v[i] + (t1 + 2.0*t2 + 2.0*t3 + t4) / 6.0;
		}
		//
		//
		bbp = y[n] - 0.0;
		ccp = w[n] - 0.0;
		eep = v[n] - 0.0;
        //
        //
        u[0] = 0.0;
		y[0] = 0.0;
		z[0] = eta;
		w[0] = 1.0;
		p[0] = ep;
		m[0] = 1.0;
		if (fabs(ha) <= 1.0) hv = sqrt(EPS);
		if (fabs(ha) > 1.0) hv = ha*sqrt(EPS);
		v[0] = ha + hv;
		    //
		//
			for (i = 0; i<n; i++)

		{
			a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			r1 = f6(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			t1 = f7(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
            //
			a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			r2 = f6(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			t2 = f7(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			//
			//		f1 returned	0.25000000000000000	double
			//
			a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			r3 = f6(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			t3 = f7(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			//
			//
		
			a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			r4 = f6(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			t4 = f7(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			//
			u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
			y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
			z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
			w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
			p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;
			m[i + 1] = m[i] + (r1 + 2.0*r2 + 2.0*r3 + r4) / 6.0;
			v[i + 1] = v[i] + (t1 + 2.0*t2 + 2.0*t3 + t4) / 6.0;
		}
		//
		//
		bbv = y[n] - 0.0;
		ccv = w[n] - 0.0;
		eev = v[n] - 0.0;
		//
		//
		dz = (cc*(bbp - bb) / hp - bb*(ccp - cc) / hp) / ((bbz - bb) / hz*(ccp - cc) / hp - (bbp - bb) / hp*(ccz - cc) / hz);
		dp = (cc*(bbz - bb) / hz - bb*(ccz - cc) / hz) / ((ccz - cc) / hz*(bbp - bb) / hp - (bbz - bb) / hz*(ccp - cc) / hp);
		dv = (cc*(bbz - bb) / hz - bb*(ccz - cc) / hz) / ((ccz - cc) / hz*(bbv - bb) / hv - (bbz - bb) / hz*(ccv - cc) / hv);
		//
		//
		eta = eta + dz;
		ep = ep + dp;
		ha = ha + dv;
		//
		//
		dd = fabs(bb) + fabs(cc) + fabs(ee);
		//
		scientificwork2 << "iter = " << iter << " dd= " << dd << endl;
		//
	} while (dd>EPS);
	//
	//
	//
	u[0] = 0.0;
	y[0] = 0.0;
	z[0] = eta;
	w[0] = 1.0;
	p[0] = ep;
	m[0] = 1.0;
	v[0] = ha;
	//
	//
				for (i = 0; i<n; i++)

		{
			a1 = f1(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			b1 = f2(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			c1 = f3(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			d1 = f4(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			g1 = f5(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			r1 = f6(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
			t1 = f7(x[i], u[i], y[i], z[i], w[i], p[i], m[i], v[i])*h;
            //
			a2 = f1(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			b2 = f2(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			c2 = f3(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			d2 = f4(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			g2 = f5(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			r2 = f6(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			t2 = f7(x[i] + 0.5*h, u[i] + 0.5*a1, y[i] + 0.5*b1, z[i] + 0.5*c1, w[i] + 0.5*d1, p[i] + 0.5*g1, m[i] + 0.5*r1, v[i] + 0.5*t1)*h;
			//
			//		f1 returned	0.25000000000000000	double
			//
			a3 = f1(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			b3 = f2(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			c3 = f3(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			d3 = f4(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			g3 = f5(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			r3 = f6(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			t3 = f7(x[i] + 0.5*h, u[i] + 0.5*a2, y[i] + 0.5*b2, z[i] + 0.5*c2, w[i] + 0.5*d2, p[i] + 0.5*g2, m[i] + 0.5*r2, v[i] + 0.5*t2)*h;
			//
			//
			a4 = f1(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			b4 = f2(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			c4 = f3(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			d4 = f4(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			g4 = f5(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			r4 = f6(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;
			t4 = f7(x[i] + h, u[i] + a3, y[i] + b3, z[i] + c3, w[i] + d3, p[i] + g3, m[i] + r3, v[i] + t3)*h;

			//
			u[i + 1] = u[i] + (a1 + 2.0*a2 + 2.0*a3 + a4) / 6.0;
			y[i + 1] = y[i] + (b1 + 2.0*b2 + 2.0*b3 + b4) / 6.0;
			z[i + 1] = z[i] + (c1 + 2.0*c2 + 2.0*c3 + c4) / 6.0;
			w[i + 1] = w[i] + (d1 + 2.0*d2 + 2.0*d3 + d4) / 6.0;
			p[i + 1] = p[i] + (g1 + 2.0*g2 + 2.0*g3 + g4) / 6.0;
			m[i + 1] = m[i] + (r1 + 2.0*r2 + 2.0*r3 + r4) / 6.0;
			v[i + 1] = v[i] + (t1 + 2.0*t2 + 2.0*t3 + t4) / 6.0;
		}
	//
	//
	for (i = 0; i<n + 1; i++)
	{
		scientificwork2<<"x[i] = "<<x[i]<<"then u[i]= "<<u[i]<<endl;
		scientificwork2 << "x[i] = " << x[i] << "then y[i]= " << y[i] << endl;
		scientificwork2 << "x[i] = " << x[i] << "then z[i]= " << z[i] << endl;
		scientificwork2 << "x[i] = " << x[i] << "then w[i]= " << w[i] << endl;
		scientificwork2 << "x[i] = " << x[i] << "then p[i]= " << p[i] << endl;
		scientificwork2 << "x[i] = " << x[i] << "then m[i]= " << m[i] << endl;
		scientificwork2 << "x[i] = " << x[i] << "then v[i]= " << v[i] << endl;

	}
	delete[]x;
	delete[]u;
	delete[]y;
	delete[]z;
	delete[]w;
	delete[]p;
	delete[]m;
	delete[]v;
	//
	scientificwork2.close();
	system("pause");
}



Thanks
Best

Is This A Good Question/Topic? 0
  • +

Page 1 of 1