More subroutines
for Stella Models on the Web

Anne Thissen
The Shodor Education Foundation
June 1998

# function of y and delta t

sub euler {
		local($y, $dt, $othervars)=@_;
		$y=$y+&function($y, $dt, $othervars);
		return($y);
}
sub heun {
		local($y, $dt, $othervars)=@_;
		$y_a=&function($y, $dt, $othervars);
		$y_b=&function($y+$y_a, $dt, $othervars);
		$y=$y+.5*$y_a+.5*$y_b;
		return($y);
}
sub runge_kutta_2 {
		local($y, $dt, $othervars)=@_;
		$y_a=&function($y, .5*$dt, $othervars);
		$y_b=&function($y+$y_a, .5*$dt, $othervars);
		$y=$y+$y_a+$y_b;
		return($y);
}
sub runge_kutta_4 {
		local($y, $dt, $othervars)=@_;
		$y_a=&function($y, $dt, $othervars);
		$y_b=&function($y+$y_a, $dt, $othervars);
		$y_c=&function($y+$y_b, $dt, $othervars);
		$y_d=&function($y+$y_c, $dt, $othervars);
		$y=$y+($y_a+2*$y_b+2*$y_c+$y_d)/6;
		return($y);
}


# y=f(y,z,dt) z=g(y,z,dt)

sub euler {
		local($y, $z, $dt, $othervars)=@_;
		$y=$y+&function($y, $z, $dt, $othervars);
		$z=$y+&function2($y, $z, $dt, $othervars);
		return($y, $z);
}
sub heun {
		local($y, $z, $dt, $othervars)=@_;
		$y_a=&function($y, $z, $dt, $othervars);
		$z_a=&function2($y+$y_a, $z, $dt, $othervars);
		$y_b=&function($y+$y_a, $z+$z_a, $dt, $othervars);
		$z_b=&function($y+$y_b, $z+$z_a, $dt, $othervars);
		$y=$y+.5*$y_a+.5*$y_b;
		$z=$z+.5*$z_a+.5*$z_b;
		return($y, $z);
}
sub runge_kutta_2 {
		local($y, $z, $dt, $othervars)=@_;
		$y_a=&function($y, $z, .5*$dt, $othervars);
		$z_a=&function2($y+$y_a, $z, .5*$dt, $othervars);
		$y_b=&function($y+$y_a, $z+$z_a, .5*$dt, $othervars);
		$z_b=&function2($y+$y_b, $z+$z_a, .5*$dt, $othervars);
		$y=$y+$y_a+$y_b;
		$z=$z+$z_a+$z_b;
		return($y, $z);
}
sub runge_kutta_4 {
		local($y, $z, $dt, $othervars)=@_;
		$y_a=&function($y, $z, $dt, $othervars);
		$z_a=&function2($y+$y_a, $z, $dt, $othervars);
		$y_b=&function($y+$y_a, $z+$z_a, $dt, $othervars);
		$z_b=&function2($y+$y_b, $z+$z_a, $dt, $othervars);
		$y_c=&function($y+$y_b, $z+$z_b, $dt, $othervars);
		$z_c=&function2($y+$y_c, $z+$z_b, $dt, $othervars);
		$y_d=&function($y+$y_c, $z+$z_c, $dt, $othervars);
		$z_d=&function2($y+$y_d, $z+$z_c, $dt, $othervars);
		$y=$y+($y_a+2*$y_b+2*$y_c+$y_d)/6;
		$z=$z+($z_a+2*$z_b+2*$z_c+$z_d)/6;
		return($y, $z);
}


# f(y,t,dt)

sub euler {
		local($y, $t, $dt, $othervars)=@_;
		$y=$y+&function($y, $t, $dt, $othervars);
		return($y);
}
sub heun {
		local($y, $t, $dt, $othervars)=@_;
		$y_a=&function($y, $t, $dt, $othervars);
		$y_b=&function($y+$y_a, $t+$dt, $dt, $othervars);
		$y=$y+.5*$y_a+.5*$y_b;
		return($y);
}
sub runge_kutta_2 {
		local($y, $t, $dt, $othervars)=@_;
		$y_a=&function($y, $t, .5*$dt, $othervars);
		$y_b=&function($y+$y_a, $t+$dt, .5*$dt, $othervars);
		$y=$y+$y_a+$y_b;
		return($y);
}
sub runge_kutta_4 {
		local($y, $t, $dt, $othervars)=@_;
		$y_a=&function($y, $t, $dt, $othervars);
		$y_b=&function($y+$y_a, $t+.5*$dt, $dt, $othervars);
		$y_c=&function($y+$y_b, $t+.5*$dt, $dt, $othervars);
		$y_d=&function($y+$y_c, $t+$dt, $dt, $othervars);
		$y=$y+($y_a+2*$y_b+2*$y_c+$y_d)/6;
		return($y);
}


# y=f(y,z,t,dt), z=g(y,z,t,dt)

sub euler {
		local($y, $z, $t, $dt, $othervars)=@_;
		$y=$y+&function($y, $z, $t, $dt, $othervars);
		$z=$z+&function2($y, $z, $t, $dt, $othervars);
		return($y, $z);
}
sub heun {
		local($y, $z, $t, $dt, $othervars)=@_;
		$y_a=&function($y, $z, $t, $dt, $othervars);
		$z_a=&function2($y+$y_a, $z, $t, $dt, $othervars);
		$y_b=&function($y+$y_a, $z+$z_a, $t+$dt, $dt, $othervars);
		$z_b=&function2($y+$y_b, $z+$z_a, $t+$dt, $dt, $othervars);
		$y=$y+.5*$y_a+.5*$y_b;
		$z=$z+.5*$z_a+.5*$z_b;
		return($y, $z);
}
sub runge_kutta_2 {
		local($y, $z, $t, $dt, $othervars)=@_;
		$y_a=&function($y, $z, $t, .5*$dt, $othervars);
		$z_a=&function2($y+$y_a, $z, $t, .5*$dt, $othervars);
		$y_b=&function($y+$y_a, $z+$z_a, $t+$dt, .5*$dt, $othervars);
		$z_b=&function2($y+$y_b, $z+$z_a, $t+$dt, .5*$dt, $othervars);
		$y=$y+$y_a+$y_b;
		return($y, $z);
}
sub runge_kutta_4 {
		local($y, $z, $t, $dt, $othervars)=@_;
		$y_a=&function($y, $z, $t, $dt, $othervars);
		$z_a=&function2($y+$y_a, $z, $t, $dt, $othervars);
		$y_b=&function($y+$y_a, $z+$z_a, $t+.5*$dt, $dt, $othervars);
		$z_b=&function2($y+$y_b, $z+$z_a, $t+.5*$dt, $dt, $othervars);
		$y_c=&function($y+$y_b, $z+$z_b, $t+.5*$dt, $dt, $othervars);
		$z_c=&function2($y+$y_c, $z+$z_b, $t+.5*$dt, $dt, $othervars);
		$y_d=&function($y+$y_c, $z+$z_c, $t+$dt, $dt, $othervars);
		$z_d=&function2($y+$y_d, $z+$z_c, $t+$dt, $dt, $othervars);
		$y=$y+($y_a+2*$y_b+2*$y_c+$y_d)/6;
		$z=$z+($z_a+2*$z_b+2*$z_c+$z_d)/6;
		return($y, $z);
}