ordinarydifferentialequation.java

数值分析算法源码(java) 这个学期一边学习java一边学习数值分析,因此用java写了一个数值分析算法的软件包numericalAnalysis. [说明] 适合使用者:会java的

package numericalAnalysis.ode;

import java.text.DecimalFormat;
import numericalAnalysis.function.Function;

/**
 * 常微分方程初值问题数值解法.对照书:数值分析第二版,史万明等编,北京理工大学出版社
 * 
 * @author 山
 * 
 */
public class OrdinaryDifferentialEquation {
	private Function function;// y'=f(x,y)的f(x,y)函数实现的接口的对象
	private double x0, y0;// 初始条件y(x0)=y0
	private DecimalFormat fmt = new DecimalFormat("0.######");

	/**
	 * 新建一个常微分方程初值问题:y'=f(x,y),x属于某区间,y(x0)=y0.
	 * 
	 * @param function
	 *            y'=f(x,y)的f(x,y)函数实现的接口的对象
	 * @param x0
	 *            初始条件y(x0)=y0
	 * @param y0
	 *            初始条件y(x0)=y0
	 */
	public OrdinaryDifferentialEquation(Function function, double x0, double y0) {
		this.function = function;
		this.x0 = x0;
		this.y0 = y0;
	}

	/**
	 * 求数值解xi,yi(i=0,1,...,n)(用龙格库塔法)
	 * 
	 * @param h
	 *            步长,则构造节点xi=x0+i*h
	 * @param xLimit
	 *            节点xi最大值(xi上限)
	 * @return 一个二维数组,其第一维索引从0到n(n为节点个数减1),代表一个节点的xi和yi对;
	 *         第二维索引只有0,1,0表示xi,1表示yi.例如:若访问这个二维数组的[3][0]索引,则访问x3,
	 *         若访问这个二维数组的[3][1]索引,则访问y3
	 */
	public double[][] getSolution(double h, double xLimit) {
		double k1, k2, k3, k4;// 对应于公式中的k1,k2,k3,k4
		double[] x, y;// 节点数组
		int n;// 节点个数减1

		n = (int) ((xLimit - x0+0.000001) / h);

		x = new double[n + 1];
		y = new double[n + 1];

		x[0] = x0;
		y[0] = y0;

		for (int i = 1; i < n + 1; i++)
			x[i] = x[0] + i * h;

		// 计算
		for (int i = 1; i < n + 1; i++) {
			k1 = h * function.f(x[i - 1], y[i - 1]);
			k2 = h * function.f(x[i - 1] + h / 2, y[i - 1] + k1 / 2);
			k3 = h * function.f(x[i - 1] + h / 2, y[i - 1] + k2 / 2);
			k4 = h * function.f(x[i - 1] + h, y[i - 1] + k3);
			y[i] = y[i - 1] + 1.0 / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
		}

		double[][] solution = new double[n + 1][2];

		for (int i = 0; i < n + 1; i++) {
			solution[i][0] = x[i];
			solution[i][1] = y[i];
		}

		return solution;
	}

	/**
	 * 欧拉公式
	 * 
	 * @param h
	 *            步长,则构造节点xi=x0+i*h
	 * @param xLimit
	 *            节点xi最大值(xi上限)
	 * @return 结果
	 */
	public String euler(double h, double xLimit){
		String output = "\n";
		double[] x, y;// 节点数组
		int n;// 节点个数减1
		
		n = (int) ((xLimit - x0+0.000001) / h);

		x = new double[n + 1];
		y = new double[n + 1];
		x[0] = x0;
		y[0] = y0;
		for (int i = 1; i < n + 1; i++)
			x[i] = x[0] + i * h;
		
		// 计算
		for (int i = 1; i < n + 1; i++) {
			y[i] = y[i - 1] + h*function.f(x[i-1],y[i-1]);
		}
		
		// 输出
		output += "xn\tyn\n";
		for (int i = 0; i < n + 1; i++)
			output += fmt.format(x[i]) + "\t" + fmt.format(y[i]) + "\n";

		return output;
	}
	
	/**
	 * 四阶古典龙格库塔法,求数值解
	 * 
	 * @param h
	 *            步长,则构造节点xi=x0+i*h
	 * @param xLimit
	 *            节点xi最大值(xi上限)
	 * @return 结果
	 */
	public String rungeKutta(double h, double xLimit) {
		String output = "\n";
		double k1, k2, k3, k4;// 对应于公式中的k1,k2,k3,k4
		double[] x, y;// 节点数组
		int n;// 节点个数减1

		n = (int) ((xLimit - x0+0.000001) / h);

		x = new double[n + 1];
		y = new double[n + 1];

		x[0] = x0;
		y[0] = y0;

		for (int i = 1; i < n + 1; i++)
			x[i] = x[0] + i * h;

		// 计算
		for (int i = 1; i < n + 1; i++) {
			k1 = h * function.f(x[i - 1], y[i - 1]);
			k2 = h * function.f(x[i - 1] + h / 2, y[i - 1] + k1 / 2);
			k3 = h * function.f(x[i - 1] + h / 2, y[i - 1] + k2 / 2);
			k4 = h * function.f(x[i - 1] + h, y[i - 1] + k3);
			y[i] = y[i - 1] + 1.0 / 6 * (k1 + 2 * k2 + 2 * k3 + k4);
		}

		// 输出
		output += "xn\tyn\n";
		for (int i = 0; i < n + 1; i++)
			output += fmt.format(x[i]) + "\t" + fmt.format(y[i]) + "\n";

		return output;
	}

	/**
	 * 显式四阶阿当姆斯公式(用龙格库塔法求表头y1,y,2,y3),求数值解
	 * 
	 * @param h
	 *            步长,则构造节点xi=x0+i*h
	 * @param xLimit
	 *            节点xi最大值(xi上限)
	 * @return 结果
	 */
	public String explicitAdams(double h, double xLimit) {
		String output = "\n";
		double[] x, y;// 节点数组
		int n;// 节点个数减1

		output+="s";
		n = (int) ((xLimit - x0+0.000001) / h);

		x = new double[n + 1];
		y = new double[n + 1];

		x[0] = x0;
		y[0] = y0;

		for (int i = 1; i < n + 1; i++)
			x[i] = x[0] + i * h;

		double[][] getHead = getSolution(h, x[3]);// 先获取表头y1,y2,y3
		for (int i = 1; i < 4; i++)
			y[i] = getHead[i][1];

		for (int i = 4; i < n + 1; i++)
			y[i] = y[i - 1]
					+ h
					/ 24
					* (55 * function.f(x[i - 1], y[i - 1]) - 59
							* function.f(x[i - 2], y[i - 2]) + 37
							* function.f(x[i - 3], y[i - 3]) - 9 * function.f(
							x[i - 4], y[i - 4]));
		
		// 输出
		output+="注意! 表头y1,y2,y3使用龙格库塔法计算,y4开始使用显式四阶阿当姆斯公式\n";
		output += "xn\tyn\n";
		for (int i = 0; i < n + 1; i++)
			output += fmt.format(x[i]) + "\t" + fmt.format(y[i]) + "\n";
		
		return output;
	}

	/**
	 * 四阶阿当姆斯预测-校正法(用龙格库塔法求表头y1,y,2,y3),求数值解
	 * 
	 * @param h
	 *            步长,则构造节点xi=x0+i*h
	 * @param xLimit
	 *            节点xi最大值(xi上限)
	 * @return 结果
	 */
	public String predictorCorrectorAdams(double h, double xLimit){
		String output = "\n";
		double[] x, y;// 节点数组
		int n;// 节点个数减1

		output+="s";
		n = (int) ((xLimit - x0+0.000001) / h);

		x = new double[n + 1];
		y = new double[n + 1];

		x[0] = x0;
		y[0] = y0;

		for (int i = 1; i < n + 1; i++)
			x[i] = x[0] + i * h;

		double[][] getHead = getSolution(h, x[3]);// 先获取表头y1,y2,y3
		for (int i = 1; i < 4; i++)
			y[i] = getHead[i][1];

		//求y4,y5,...,yn
		for (int i = 3; i < n; i++){
			//求y[i+1]*
			double yStar;
			yStar = y[i]
					+ h
					/ 24
					* (55 * function.f(x[i], y[i]) - 59
							* function.f(x[i - 1], y[i - 1]) + 37
							* function.f(x[i - 2], y[i - 2]) - 9 * function.f(
							x[i - 3], y[i - 3]));
			
			y[i+1] = y[i]
						+ h
						/ 24
						* (9 * function.f(x[i + 1], yStar) + 19
								* function.f(x[i], y[i]) -5
								* function.f(x[i - 1], y[i - 1]) + function.f(
								x[i - 2], y[i - 2]));
		}
		
		// 输出
		output+="注意! 表头y1,y2,y3使用龙格库塔法计算,y4开始使用显式四阶阿当姆斯公式\n";
		output += "xn\tyn\n";
		for (int i = 0; i < n + 1; i++)
			output += fmt.format(x[i]) + "\t" + fmt.format(y[i]) + "\n";
		
		return output;
	}
	
}