//  ------------------------------------------------------
//  -------  Prüfung 4.9      Ib03        ----------------
//  -------                               ----------------
//  -------  Runge-Kuttta für DG-System   ----------------
//  -------  2. Ordnung                   ----------------
//  -------                               ----------------
//  -------  21.04.2005    M. Vogel       ----------------
//  ------------------------------------------------------

public class Runge_Kutta_System
{
	public static double y10 = 4700;	       // Anfangsbedingung 1. Gleichung
	public static double y20 = 3000;	       // Anfangsbedingung 2. Gleichung
	public static double t0  = 0;
	public static double tEnd= 10;
	public static int    n   = 100;

	// Arrays
	public static double []   t = new double[n+1];
	public static double [][] w = new double[3][n+1];	// 1. für w1, 2. für w2
	
	// Funktionen
	public static double f1(double t, double y1, double y2) { return y1 * (4 - 0.0003*y1 - 0.0004*y2) ; }
	public static double f2(double t, double y1, double y2) { return y2 * (2 - 0.0002*y1 - 0.0001*y2) ; }
	
	public static void main (String [] args)
	{
		System.out.println((tEnd - t0)/n);
		//getInfo();
		System.out.println("\n\n");
		proc(y10, y20, t0, tEnd, n);
		print();
	}
	
	public static void proc(double y10, double y20, double t0, double tEnd, int n)
	{
		//double h = (tEnd - t0)/n;
		double h = 0.1;
		t[0]     = t0;
		w[1][0]  = y10;
		w[2][0]  = y20;
		
		for(int j=0; j<n; j++)
		{
			double k11 = h * f1( t[j]           , w[1][j]           , w[2][j]);
			double k12 = h * f2( t[j]           , w[1][j]           , w[2][j]);
			double k21 = h * f1( t[j] + h/2.0   , w[1][j]+k11/2.0   , w[2][j]+k12/2.0);
			double k22 = h * f2( t[j] + h/2.0   , w[1][j]+k11/2.0   , w[2][j]+k12/2.0);
			double k31 = h * f1( t[j] + h/2.0   , w[1][j]+k21/2.0   , w[2][j]+k22/2.0);
			double k32 = h * f2( t[j] + h/2.0   , w[1][j]+k21/2.0   , w[2][j]+k22/2.0);
			double k41 = h * f1( t[j]+h         , w[1][j]+k31       , w[2][j]+k32); 
			double k42 = h * f2( t[j]+h         , w[1][j]+k31       , w[2][j]+k32); 
			t[j+1] = t[j]+h; 
			w[1][j+1]  = w[1][j] + 1.0/6.0*(k11 + 2*k21 +2*k31 + k41);
			w[2][j+1]  = w[2][j] + 1.0/6.0*(k12 + 2*k22 +2*k32 + k42);
		}
	}
	
	public static void print(){
		System.out.println("Ausgabe: \n");
		for(int i = 0; i<w[0].length; i++)
		{
			System.out.print("x: ");
			InOut.print(t[i], 2,2);
			//System.out.println("x: " + t[i] +'\t'+'\t' +"w[1]["+i+"]: " + w[1][i] +'\t'+'\t'+"w[2]["+i+"]: "+ w[2][i]);
			System.out.println("\t\tw[1]["+i+"]: " + w[1][i] +'\t'+'\t'+"w[2]["+i+"]: "+ w[2][i]);
		}
	}
}