Package to solve ODE with Newton Raphson and nested Newton algorithm

Author: niccolotubini

src/main/java/it/geoframe/blogspot/numerical/ode/NestedNewton.java

@@ -0,0 +1,83 @@
+/*
+ * GNU GPL v3 License
+ *
+ * Copyright 2021 Niccolò Tubini, Giuseppe Formetta, Riccardo Rigon
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+package it.geoframe.blogspot.numerical.ode;
+
+
+/**
+ * @author Niccolò Tubini, Giuseppe Formetta
+ * 
+ */
+public class NestedNewton {
+
+	private double x0;
+	
+	private double x_k;
+	
+	private double tol = 1e-10;
+	
+	private double error;
+	
+	private int iteration;
+	
+	private int maxIteration = 100;
+	
+	private double f;
+	
+	private double ff;
+	
+	public double solve( double x, OrdinaryDifferentialEquation ode ) {
+		
+		this.x0 = x;
+		
+		error = 1.0;
+		
+		iteration = 1;
+				
+		
+		for(int outerIteration=0; outerIteration<maxIteration; outerIteration++) {
+			
+			f = ode.compute(x0) - ode.computeRHS();
+			
+			if(Math.abs(f)<tol) {
+				
+				return x0;
+				
+			}
+			
+			x_k = x0;
+			
+			for(int innerrIteration=0; innerrIteration<maxIteration; innerrIteration++) {
+				
+				ff = ode.computePIntegral(x0) - (ode.computeQIntegral(x_k)+ode.computeQ(x_k)*(x0-x_k)) - ode.computeRHS();
+				
+				if(Math.abs(ff)<tol) {
+					
+					return x0;
+					
+				}
+				
+				x0 = x0 - 1/(ode.computeP(x0)-ode.computeQ(x_k))*ff;
+			}
+		}
+
+		return x0;
+	}
+	
+}

src/main/java/it/geoframe/blogspot/numerical/ode/NewtonRaphson.java

@@ -0,0 +1,74 @@
+/*
+ * GNU GPL v3 License
+ *
+ * Copyright 2021 Niccolò Tubini, Giuseppe Formetta, Riccardo Rigon
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+package it.geoframe.blogspot.numerical.ode;
+
+
+/**
+ * @author Niccolò Tubini, Giuseppe Formetta
+ * 
+ */
+public class NewtonRaphson {
+
+	private double x0;
+	
+	private double xNew;
+	
+	private double tol = 1e-10;
+	
+	private double error;
+	
+	private int iteration;
+	
+	private int maxIteration = 100;
+	
+	
+	
+	public double solve( double x, OrdinaryDifferentialEquation ode ) {
+		
+		this.x0 = x;
+		
+		error = 1.0;
+		
+		iteration = 1;
+		
+		
+		while(error>tol) {
+		
+			if(iteration>maxIteration) {
+				
+				System.out.println("\t NewtonRaphson reach maximum numeber of iteration " + ode.toString());
+				break;
+				
+			}
+			
+			xNew = -ode.compute(this.x0)/ode.computeDerivative(this.x0) + this.x0;
+			
+			error = Math.abs(xNew-this.x0)/this.x0;
+			
+			this.x0 = xNew;
+			
+			iteration++;
+				
+		}
+		
+		return xNew;
+	}
+	
+}

src/main/java/it/geoframe/blogspot/numerical/ode/ODELiquidMass.java

@@ -0,0 +1,87 @@
+/*
+ * GNU GPL v3 License
+ *
+ * Copyright 2021 Niccolò Tubini, Giuseppe Formetta, Riccardo Rigon
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+package it.geoframe.blogspot.numerical.ode;
+
+/**
+ * @author Niccolò Tubini, Giuseppe Formetta
+ * 
+ */
+public class ODELiquidMass implements OrdinaryDifferentialEquation {
+
+	private double initialCondition;
+	private double rainfall;
+	private double freezing;
+	private double melting;
+	private double snowPorosity;
+	
+	
+	public ODELiquidMass(double initialCondition, double rainfall, double freezing, double melting, double snowPorosity) {
+		
+		this.initialCondition = initialCondition;
+		this.rainfall = rainfall;
+		this.freezing = freezing;
+		this.melting = melting;
+		this.snowPorosity = snowPorosity;
+		
+	}
+	
+	@Override
+	public double compute(double liquidWater) {
+		// TODO Auto-generated method stub
+
+		return Math.min(liquidWater,snowPorosity) + Math.max(liquidWater-snowPorosity, 0) - initialCondition - rainfall + freezing - melting;
+		
+	}
+
+	@Override
+	public double computeDerivative(double liquidWater) {
+		// TODO Auto-generated method stub
+		
+//		if(liquidWater>snowPorosity) {
+//			
+//			return 2;
+//			
+//		} else {
+			
+			return 1;
+			
+//		}
+		
+	}
+
+	@Override
+	public double computeP(double x) {
+		// TODO Auto-generated method stub
+		return 0;
+	}
+
+	@Override
+	public double computePIntegral(double x) {
+		// TODO Auto-generated method stub
+		return 0;
+	}
+
+	@Override
+	public double computeRHS() {
+		// TODO Auto-generated method stub
+		return 0;
+	}
+
+}

src/main/java/it/geoframe/blogspot/numerical/ode/ODESolidMass.java

@@ -0,0 +1,74 @@
+/*
+ * GNU GPL v3 License
+ *
+ * Copyright 2021 Niccolò Tubini, Giuseppe Formetta, Riccardo Rigon
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+package it.geoframe.blogspot.numerical.ode;
+
+/**
+ * @author Niccolò Tubini, Giuseppe Formetta
+ * 
+ */
+public class ODESolidMass implements OrdinaryDifferentialEquation {
+
+	private double initialCondition;
+	private double snowfall;
+	private double freezing;
+	private double melting;
+	
+	
+	public ODESolidMass(double initialCondition, double snowfall, double freezing, double melting) {
+		
+		this.initialCondition = initialCondition;
+		this.snowfall = snowfall;
+		this.freezing = freezing;
+		this.melting = melting;
+		
+	}
+	
+	@Override
+	public double compute(double solidWater) {
+		// TODO Auto-generated method stub
+		
+		return solidWater - initialCondition - snowfall - freezing + melting;
+	}
+
+	@Override
+	public double computeDerivative(double x) {
+		// TODO Auto-generated method stub
+		return 1;
+	}
+
+	@Override
+	public double computeP(double x) {
+		// TODO Auto-generated method stub
+		return 0;
+	}
+
+	@Override
+	public double computePIntegral(double x) {
+		// TODO Auto-generated method stub
+		return 0;
+	}
+
+	@Override
+	public double computeRHS() {
+		// TODO Auto-generated method stub
+		return 0;
+	}
+
+}

src/main/java/it/geoframe/blogspot/numerical/ode/OrdinaryDifferentialEquation.java

@@ -0,0 +1,52 @@
+/*
+ * GNU GPL v3 License
+ *
+ * Copyright 2021 Niccolò Tubini, Giuseppe Formetta, Riccardo Rigon
+ *
+ * This program is free software: you can redistribute it and/or modify
+ * it under the terms of the GNU General Public License as published by
+ * the Free Software Foundation, either version 3 of the License, or
+ * (at your option) any later version.
+ *
+ * This program is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ * GNU General Public License for more details.
+ *
+ * You should have received a copy of the GNU General Public License
+ * along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+package it.geoframe.blogspot.numerical.ode;
+
+/**
+ * @author Niccolò Tubini, Giuseppe Formetta
+ * 
+ */
+
+public interface OrdinaryDifferentialEquation {
+	
+		
+	public abstract double compute(double x);
+	
+	public abstract double computeDerivative(double x);
+	
+	public abstract double computeP(double x);
+		
+	public abstract double computePIntegral(double x);
+	
+	public abstract double computeRHS();
+	
+	public default  double computeQ(double x) {
+		
+		return computeP(x) - computeDerivative(x);
+				
+	}
+
+	public default double computeQIntegral(double x) {
+		
+		return computePIntegral(x) - compute(x);
+		
+	}
+
+}