Here is a listing of the most important methods used in the
JBONE applet; to obtain a complete version of the
source code, qualifying users are kindly requested to use
the download service.
jbone: /* $Id: jbone.java,v 1.46 2001/01/16 18:35:55 andrej Exp $ */
jbone:
jbone: import java.util.*;
jbone: import java.text.*;
jbone: import java.awt.*;
jbone: import java.lang.Math;
jbone: import java.applet.Applet;
jbone:
jbone:
jbone:
jbone: /*****************************************************************************
jbone: * Java Bed for ONE dimensional partial differential equations.
jbone: *
jbone: * Andre JAUN (jaun@kth.se) and Johan HEDIN (johanh@fusion.kth.se)
jbone: * Royal Institute of Technology, SE-100 44 Stockholm.
jbone: * Copyright 1998-2004. All Rights Reserved.
jbone: * @version 4.2 (CVS/RCS $Revision: 1.47 $)
jbone: * Educational software distributed by www.lifelong-learners.com.
jbone: *
jbone: * This shareware can be obtained by e-mail from the authors.
jbone: * Permission to use, copy, and modify the source and its documentation
jbone: * for your own personal use is hereby granted provided that this copyright
jbone: * notice appears in all copies. The authors makes no representations or
jbone: * warranties about the suitability of the software, either express or
jbone: * implied, including but not limited to the implied warranties of
jbone: * merchantability, fitness for a particular purpose, or non-infringement.
jbone: * The authors shall not be liable for any damages suffered by licensee
jbone: * as a result of using or modifying this software or its derivatives.
jbone: ******************************************************************************/
jbone: public class jbone extends Applet implements RunDataNotable, Runnable {
jbone:
jbone: // Constant for the different PDEs, inital condition, methods and schemes
jbone:
jbone:
jbone: /** The advection diffusion equation */
jbone: final static String ADVECTION = "Advection";
jbone: /** The wave equation */
jbone: final static String WAVE = "Wave";
jbone: /** The Burger equation for shock waves */
jbone: final static String BURGER = "Burger (shock)";
jbone: /** The Korteweg-De Vrie's equation for solitons */
jbone: final static String KDV = "KdV (solitons)";
jbone: /** The Schroedinger equation in quantum mechanics */
jbone: final static String SCHROED = "Schroedinger";
jbone: /** The Black-Scholes equation for option pricing */
jbone: final static String BSCHOLES = "Black-Scholes";
jbone: /** Student equation 1 */
jbone: final static String MY_PDE_1 = "My project 1";
jbone: /** Student equation 2 */
jbone: final static String MY_PDE_2 = "My project 2";
jbone: /** Solution of the exercises */
jbone: final static String EXERCISE = "Exercise";
jbone: /** Vector of the names of the PDE types @see Solution */
jbone: final static String[] pdeNames = {ADVECTION, WAVE, BURGER, KDV, SCHROED,
jbone: BSCHOLES, MY_PDE_1, MY_PDE_2, EXERCISE };
jbone: /** GUI list of all the pdes */
jbone: private MyChoice pdeSelection;
jbone:
jbone:
jbone: /** Solve with finite differances */
jbone: final static String FD = "Finite differences";
jbone: /** Solve with finite elements method */
jbone: final static String FEM = "Finite elements";
jbone: /** Solve with Fourier transformations */
jbone: final static String FFT = "Fourier transform";
jbone: /** Solve with Lagrangian scheme */
jbone: final static String CHA = "Lagrangian methods";
jbone: /** Solve with particle methods */
jbone: final static String PPP = "Particle methods";
jbone: /** Solve with particle methods */
jbone: final static String PPPS= "Particle methods*";
jbone: /** Solve with sampling methods */
jbone: final static String SMP = "Sampling methods";
jbone: /** Solve with sampling methods */
jbone: final static String SMPS= "Sampling methods*";
jbone: /** Vector of the names of all the solvers @see Solution */
jbone: static String[] methodNames = {FD, FEM, FFT, CHA, PPP, PPPS, SMP, SMPS};
jbone: /** GUI list of all the methods */
jbone: private MyChoice methodSelection;
jbone:
jbone:
jbone: /** Solve with Standard scheme */
jbone: final static String DEFAULT = "Standard scheme";
jbone: /** Solve with Explicit 2-level */
jbone: final static String EXP2LVL = "Explicit 2-level";
jbone: /** Solve with Explicit 3-level */
jbone: final static String EXP3LVL = "Explicit 3-level";
jbone: /** Solve with Implicit 2-level */
jbone: final static String IMP2LVL = "Implicit 2-level";
jbone: /** Solve with Expl-LaxWendroff */
jbone: final static String EXPLAX = "Expl-LaxWendroff";
jbone: /** Solve with Impl-LaxWendroff */
jbone: final static String IMPLAX = "Impl-LaxWendroff";
jbone: /** Solve with Lax-Friedrichs */
jbone: final static String LAXFRIED = "Lax-Friedrichs";
jbone: /** Solve with Leap-frog (FDTD) */
jbone: final static String LPFROG = "Leap-frog (FDTD)";
jbone: /** FDTD with out-going and perfectly matched layer BC */
jbone: final static String FDTDABS = "FDTD <out--pml>";
jbone: /** European put without normalizing */
jbone: final static String EUNAIVE = "European naive";
jbone: /** European put in normalized variables */
jbone: final static String EUPUT = "European vanilla";
jbone: /** American put without normalizing */
jbone: final static String AMNAIVE = "American naive";
jbone: /** American put in normalized variables */
jbone: final static String AMPUT = "American vanilla";
jbone: /** Solve with Tunable Integration */
jbone: final static String TUNFEM = "Tune Integration";
jbone: /** Solve with Expanded convolution */
jbone: final static String EXPAND = "Expanded convol.";
jbone: /** Solve with Aliased convolution */
jbone: final static String ALIASED = "Aliased convol.";
jbone: /** Solve with CubicHermite FEM */
jbone: final static String CUBFEM = "CubicHermite FEM";
jbone: /** Solve with Cubic -- Splines */
jbone: final static String CUBSPL = "Cubic---Splines";
jbone: /** Solve with Monte-Carlo using a normally distributed random walkers */
jbone: final static String MCNORM = "MC norm. walk";
jbone: /** Solve with Monte-Carlo using log-normally distributed random walkers */
jbone: final static String MCLOGN = "MC log-normal";
jbone: /** Solve by sampling with a tree */
jbone: final static String TREE2 = "Binomial tree";
jbone: /** Solve with My scheme */
jbone: final static String MY_SCHEME = "My scheme";
jbone: /** Solve with My scheme 2 */
jbone: final static String MY_SCHEME_2 = "My scheme 2";
jbone: /** Solution of exercise 1.01 */ final static String EXC1_01="Exercise 1.01";
jbone: /** Solution of exercise 1.02 */ final static String EXC1_02="Exercise 1.02";
jbone: /** Solution of exercise 1.03 */ final static String EXC1_03="Exercise 1.03";
jbone: /** Solution of exercise 1.04 */ final static String EXC1_04="Exercise 1.04";
jbone: /** Solution of exercise 1.05 */ final static String EXC1_05="Exercise 1.05";
jbone: /** Solution of exercise 1.06 */ final static String EXC1_06="Exercise 1.06";
jbone: /** Solution of exercise 1.07 */ final static String EXC1_07="Exercise 1.07";
jbone: /** Solution of exercise 1.08 */ final static String EXC1_08="Exercise 1.08";
jbone: /** Solution of exercise 1.09 */ final static String EXC1_09="Exercise 1.09";
jbone: /** Solution of exercise 1.10 */ final static String EXC1_10="Exercise 1.10";
jbone: /** Solution of exercise 2.01 */ final static String EXC2_01="Exercise 2.01";
jbone: /** Solution of exercise 2.02 */ final static String EXC2_02="Exercise 2.02";
jbone: /** Solution of exercise 2.03 */ final static String EXC2_03="Exercise 2.03";
jbone: /** Solution of exercise 2.04 */ final static String EXC2_04="Exercise 2.04";
jbone: /** Solution of exercise 2.05 */ final static String EXC2_05="Exercise 2.05";
jbone: /** Solution of exercise 2.06 */ final static String EXC2_06="Exercise 2.06";
jbone: /** Solution of exercise 2.07 */ final static String EXC2_07="Exercise 2.07";
jbone: /** Solution of exercise 2.08 */ final static String EXC2_08="Exercise 2.08";
jbone: /** Solution of exercise 2.09 */ final static String EXC2_09="Exercise 2.09";
jbone: /** Solution of exercise 2.10 */ final static String EXC2_10="Exercise 2.10";
jbone: /** Solution of exercise 3.01 */ final static String EXC3_01="Exercise 3.01";
jbone: /** Solution of exercise 3.02 */ final static String EXC3_02="Exercise 3.02";
jbone: /** Solution of exercise 3.03 */ final static String EXC3_03="Exercise 3.03";
jbone: /** Solution of exercise 3.04 */ final static String EXC3_04="Exercise 3.04";
jbone: /** Solution of exercise 3.05 */ final static String EXC3_05="Exercise 3.05";
jbone: /** Solution of exercise 3.06 */ final static String EXC3_06="Exercise 3.06";
jbone: /** Solution of exercise 3.07 */ final static String EXC3_07="Exercise 3.07";
jbone: /** Solution of exercise 3.08 */ final static String EXC3_08="Exercise 3.08";
jbone: /** Solution of exercise 3.09 */ final static String EXC3_09="Exercise 3.09";
jbone: /** Solution of exercise 3.10 */ final static String EXC3_10="Exercise 3.10";
jbone: /** Solution of exercise 4.01 */ final static String EXC4_01="Exercise 4.01";
jbone: /** Solution of exercise 4.02 */ final static String EXC4_02="Exercise 4.02";
jbone: /** Solution of exercise 4.03 */ final static String EXC4_03="Exercise 4.03";
jbone: /** Solution of exercise 4.04 */ final static String EXC4_04="Exercise 4.04";
jbone: /** Solution of exercise 4.05 */ final static String EXC4_05="Exercise 4.05";
jbone: /** Solution of exercise 4.06 */ final static String EXC4_06="Exercise 4.06";
jbone: /** Solution of exercise 4.07 */ final static String EXC4_07="Exercise 4.07";
jbone: /** Solution of exercise 4.08 */ final static String EXC4_08="Exercise 4.08";
jbone: /** Solution of exercise 4.09 */ final static String EXC4_09="Exercise 4.09";
jbone: /** Solution of exercise 4.10 */ final static String EXC4_10="Exercise 4.10";
jbone: /** Solution of exercise 5.01 */ final static String EXC5_01="Exercise 5.01";
jbone: /** Solution of exercise 5.02 */ final static String EXC5_02="Exercise 5.02";
jbone: /** Solution of exercise 5.03 */ final static String EXC5_03="Exercise 5.03";
jbone: /** Solution of exercise 5.04 */ final static String EXC5_04="Exercise 5.04";
jbone: /** Solution of exercise 5.05 */ final static String EXC5_05="Exercise 5.05";
jbone: /** Solution of exercise 5.06 */ final static String EXC5_06="Exercise 5.06";
jbone: /** Solution of exercise 5.07 */ final static String EXC5_07="Exercise 5.07";
jbone: /** Solution of exercise 5.08 */ final static String EXC5_08="Exercise 5.08";
jbone: /** Solution of exercise 5.09 */ final static String EXC5_09="Exercise 5.09";
jbone: /** Solution of exercise 5.10 */ final static String EXC5_10="Exercise 5.10";
jbone: /** Solution of exercise 6.01 */ final static String EXC6_01="Exercise 6.01";
jbone: /** Solution of exercise 6.02 */ final static String EXC6_02="Exercise 6.02";
jbone: /** Solution of exercise 6.03 */ final static String EXC6_03="Exercise 6.03";
jbone: /** Solution of exercise 6.04 */ final static String EXC6_04="Exercise 6.04";
jbone: /** Solution of exercise 6.05 */ final static String EXC6_05="Exercise 6.05";
jbone: /** Solution of exercise 6.06 */ final static String EXC6_06="Exercise 6.06";
jbone: /** Solution of exercise 6.07 */ final static String EXC6_07="Exercise 6.07";
jbone: /** Solution of exercise 6.08 */ final static String EXC6_08="Exercise 6.08";
jbone: /** Solution of exercise 6.09 */ final static String EXC6_09="Exercise 6.09";
jbone: /** Solution of exercise 6.10 */ final static String EXC6_10="Exercise 6.10";
jbone: /** Solution of exercise 7.01 */ final static String EXC7_01="Exercise 7.01";
jbone: /** Solution of exercise 7.02 */ final static String EXC7_02="Exercise 7.02";
jbone: /** Solution of exercise 7.03 */ final static String EXC7_03="Exercise 7.03";
jbone: /** Solution of exercise 7.04 */ final static String EXC7_04="Exercise 7.04";
jbone: /** Solution of exercise 7.05 */ final static String EXC7_05="Exercise 7.05";
jbone: /** Solution of exercise 7.06 */ final static String EXC7_06="Exercise 7.06";
jbone: /** Solution of exercise 7.07 */ final static String EXC7_07="Exercise 7.07";
jbone: /** Solution of exercise 7.08 */ final static String EXC7_08="Exercise 7.08";
jbone: /** Solution of exercise 7.09 */ final static String EXC7_09="Exercise 7.09";
jbone: /** Solution of exercise 7.10 */ final static String EXC7_10="Exercise 7.10";
jbone: /** Solution of exercise A.01 */ final static String EXCA_01="Exercise A.01";
jbone: /** Solution of exercise B.01 */ final static String EXCB_01="Exercise B.01";
jbone: /** Solution of exercise C.01 */ final static String EXCC_01="Exercise C.01";
jbone: /** Solution of exercise D.01 */ final static String EXCD_01="Exercise D.01";
jbone: /** Solution of exercise E.01 */ final static String EXCE_01="Exercise E.01";
jbone: /** Solution of exercise F.01 */ final static String EXCF_01="Exercise F.01";
jbone: /** Vector of the names of all the schemes @see Solution */
jbone: final static String[] schemeNames =
jbone: {DEFAULT, EXP2LVL, EXP3LVL, IMP2LVL, EXPLAX, IMPLAX, LAXFRIED,
jbone: LPFROG, FDTDABS, TUNFEM, EXPAND, ALIASED, CUBFEM, CUBSPL,
jbone: EUNAIVE, AMNAIVE, EUPUT, AMPUT,
jbone: MCNORM, MCLOGN, TREE2, MY_SCHEME, MY_SCHEME_2,
jbone: EXC1_01, EXC1_02, EXC1_03, EXC1_04, EXC1_05,
jbone: EXC1_06, EXC1_07, EXC1_08, EXC1_09, EXC1_10,
jbone: EXC2_01, EXC2_02, EXC2_03, EXC2_04, EXC2_05,
jbone: EXC2_06, EXC2_07, EXC2_08, EXC2_09, EXC2_10,
jbone: EXC3_01, EXC3_02, EXC3_03, EXC3_04, EXC3_05,
jbone: EXC3_06, EXC3_07, EXC3_08, EXC3_09, EXC3_10,
jbone: EXC4_01, EXC4_02, EXC4_03, EXC4_04, EXC4_05,
jbone: EXC4_06, EXC4_07, EXC4_08, EXC4_09, EXC4_10,
jbone: EXC5_01, EXC5_02, EXC5_03, EXC5_04, EXC5_05,
jbone: EXC5_06, EXC5_07, EXC5_08, EXC5_09, EXC5_10,
jbone: EXC6_01, EXC6_02, EXC6_03, EXC6_04, EXC6_05,
jbone: EXC6_06, EXC6_07, EXC6_08, EXC6_09, EXC6_10,
jbone: EXC7_01, EXC7_02, EXC7_03, EXC7_04, EXC7_05,
jbone: EXC7_06, EXC7_07, EXC7_08, EXC7_09, EXC7_10
jbone: };
jbone: /** GUI list of all the schemes */
jbone: private MyChoice schemeSelection;
jbone:
jbone:
jbone:
jbone: /** A box as initial condition */
jbone: final static String BOX = "Box";
jbone: /** A gaussian as initial condition */
jbone: final static String GAUSSIAN = "Gaussian";
jbone: /** A cosinus as initial condition */
jbone: final static String COSINUS = "Cosine";
jbone: /** A soliton as initial condition */
jbone: final static String SOLITON = "Soliton";
jbone: /** A wave-packet as initial condition */
jbone: final static String WAVEPACKET = "WavePacket";
jbone: /** A stockmarket option as initial condition */
jbone: final static String OPTION = "PutOption";
jbone: /** Vector of the names of allowed initial conditions @see ShapeFunction */
jbone: final static String[] icNames = {BOX, GAUSSIAN, COSINUS, SOLITON,
jbone: WAVEPACKET, OPTION};
jbone: /** GUI list of all the ICs */
jbone: private MyChoice icSelection;
jbone:
jbone: /** Operate edit mode with TAG parameters displayed by default */
jbone: final static String EDIT = "Double-click below:";
jbone: /** Operate edit mode with ALL parameters displayed */
jbone: final static String EDITALL = "Show all parameters";
jbone: /** Operate a console output of the function values */
jbone: final static String CONSOLE = "Print data to console ";
jbone: /** Vector of the names of all the operations */
jbone: static String[] operNames = {EDIT, EDITALL, CONSOLE};
jbone: /** GUI list of all the operations */
jbone: private MyChoice operSelection;
jbone:
jbone: /** The plot area */
jbone: private PlotArea plotArea;
jbone: /** The run parameters */
jbone: private RunData runData;
jbone: /** Text on the Start/Stop button */
jbone: private final String startName = "Start/Stop";
jbone: /** Text on the Step 1 button */
jbone: private final String step1Name = "Step 1";
jbone: /** Text on the Step 10 button */
jbone: private final String step10Name = "Step 10";
jbone: /** Text on the toggle display button */
jbone: private final String displayName = "Toggle Display";
jbone: /** Text on the initialize button */
jbone: private final String initializeName = "Initialize";
jbone:
jbone: /** Whether the simulation is running */
jbone: private boolean frozen = true;
jbone: /** Current step number */
jbone: private int step = 0;
jbone: /** Operate nsteps before stopping */
jbone: private int nstep = -1;
jbone: /** Milliseconds between plots */
jbone: private int delay = 84;
jbone: /** Thread label */
jbone: Thread runThread;
jbone: /** Potentially reset by main */
jbone: private boolean isAnApplet = true;
jbone:
jbone: /** The solution */
jbone: private Solution solution;
jbone: /** The physical data */
jbone: private PhysData physData;
jbone:
jbone:
jbone: /** Information
jbone: ****************************************************************************/
jbone: public String getAppletInfo() {
jbone: return "JBONE 4.0 -- an educational software (C) 1997-2004 A.Jaun,J.Hedin";
jbone: }
jbone:
jbone:
jbone: /** Master initialization and layout
jbone: ****************************************************************************/
jbone: public void init() {
jbone: pdeSelection = new MyChoice(pdeNames);
jbone: methodSelection = new MyChoice(methodNames);
jbone: schemeSelection = new MyChoice(schemeNames);
jbone: icSelection = new MyChoice(icNames);
jbone: operSelection = new MyChoice(operNames);
jbone: runData = new RunData(this);
jbone:
jbone: if (isAnApplet) tagModify();
jbone: createWindow(); createPhysics();
jbone: createSolution();
jbone: setInitialCondition(solution, null, false);
jbone: } // init
jbone:
jbone: /** Instanciate a GUI window
jbone: @see Solution
jbone: ***************************************************************************/
jbone: private void createWindow(){
jbone: // First the layout
jbone: GridBagLayout gb = new GridBagLayout();
jbone: GridBagConstraints c = new GridBagConstraints();
jbone:
jbone: this.setBackground(new Color(Integer.parseInt("88FFFF",16)));
jbone: setFont(new Font("Courier", Font.PLAIN, 12));
jbone: setLayout(gb);
jbone:
jbone: // Prepare for top row of buttons
jbone: c.gridwidth = 1;
jbone: c.gridheight = 1;
jbone: c.weightx = 1;
jbone: c.fill = GridBagConstraints.HORIZONTAL;
jbone:
jbone: // Engage
jbone: gbAdd(gb, c, methodSelection);
jbone: gbAdd(gb, c, schemeSelection);
jbone: gbAdd(gb, c, icSelection);
jbone: gbAdd(gb, c, pdeSelection);
jbone: c.gridwidth = GridBagConstraints.REMAINDER; //end of row
jbone: gbAdd(gb, c, operSelection);
jbone:
jbone: // Mr Data, take us to the next row please
jbone: c.gridwidth = GridBagConstraints.RELATIVE;
jbone: c.gridheight = 1;
jbone: c.gridwidth = 4;
jbone: c.weightx = 0;
jbone: c.weighty = 1;
jbone: c.fill = GridBagConstraints.BOTH;
jbone: plotArea = new PlotArea();
jbone: gbAdd(gb, c, plotArea);
jbone: c.gridwidth = GridBagConstraints.REMAINDER; //end of row
jbone:
jbone: gbAdd(gb, c, runData);
jbone:
jbone: // Finally the last row
jbone: c.gridwidth = 1;
jbone: c.weighty = 0;
jbone: c.fill = GridBagConstraints.HORIZONTAL;
jbone: gbAdd(gb, c, new Button(startName));
jbone: gbAdd(gb, c, new Button(step1Name));
jbone: gbAdd(gb, c, new Button(step10Name));
jbone: gbAdd(gb, c, new Button(displayName));
jbone: c.gridwidth = GridBagConstraints.REMAINDER; //end of row
jbone: gbAdd(gb, c, new Button(initializeName));
jbone:
jbone: } // createWindow
jbone:
jbone:
jbone: /** Helper method for adding objects to a GridBagLayout
jbone: @param gb The layout
jbone: @param c The constraints
jbone: @param item The object to add
jbone: */
jbone: private void gbAdd(GridBagLayout gb, GridBagConstraints c, Component item){
jbone: gb.setConstraints(item, c);
jbone: add(item);
jbone: } // gbAdd
jbone:
jbone:
jbone: /** Instanciate the physical data (such as potential) from the parameters.
jbone: @see Solution
jbone: @return A new set of physical parameters
jbone: ***************************************************************************/
jbone: private void createPhysics(){
jbone: physData = new PhysData(runData);
jbone: } // createPhysics
jbone:
jbone: /** Instanciate a solution and select the method and scheme for computations
jbone: @see Solution
jbone: @return A new Solution
jbone: ***************************************************************************/
jbone: private void createSolution(){
jbone: runData.createMesh();
jbone: if(methodSelection.getSelectedItem().equals(FD)) {
jbone: solution = new FDSolution(runData);
jbone: } else if(methodSelection.getSelectedItem().equals(FEM)){
jbone: solution = new FEMSolution(runData);
jbone: } else if(methodSelection.getSelectedItem().equals(FFT)){
jbone: solution = new FFTSolution(runData);
jbone: } else if(methodSelection.getSelectedItem().equals(CHA)){
jbone: solution = new CHASolution(runData);
jbone: } else if(methodSelection.getSelectedItem().equals(PPP)){
jbone: solution = new MCPSolution(runData);
jbone: } else if(methodSelection.getSelectedItem().equals(PPPS)){
jbone: solution = new MCPDrawSolution(runData);
jbone: } else if(methodSelection.getSelectedItem().equals(SMP)){
jbone: solution = new MCSSolution(runData);
jbone: } else if(methodSelection.getSelectedItem().equals(SMPS)){
jbone: solution = new MCSDrawSolution(runData);
jbone: } else {
jbone: /** throw new MethodException("Can't find method"
jbone: + methodSelection.getSelectedItem()); */
jbone: } // if
jbone:
jbone: String ch; //Synchronize the choices is delicate
jbone: pdeSelection.sync(solution);
jbone: ch=pdeSelection.getSelectedItem();
jbone: if(ch==null) System.out.println("ERROR in jbone: pdeSelection ="+ch);
jbone: solution.setPde(ch);
jbone: ch=methodSelection.getSelectedItem();
jbone: if(ch==null) System.out.println("ERROR in jbone: methodSelection ="+ch);
jbone: solution.setMethod(ch);
jbone: schemeSelection.sync(solution);
jbone: ch=schemeSelection.getSelectedItem();
jbone: if(ch==null) System.out.println("ERROR in jbone: schemeSelection ="+ch);
jbone: solution.setScheme(ch);
jbone: plotArea.setSolution(solution);
jbone: } // createSolution
jbone:
jbone: /** Set the initial condition according to the runData parameters.
jbone: @param inSolution The solution to modify
jbone: @param oldSolution obtained previously with another method, scheme
jbone: @param keepFunction Whether to keep the old function
jbone: @see Solution
jbone: ****************************************************************************/
jbone: private void setInitialCondition(Solution inSolution,
jbone: ShapeFunction oldSolution,
jbone: boolean keepFunction){
jbone: if(keepFunction){
jbone: try{
jbone: inSolution.discretize(oldSolution);
jbone: } catch (NullPointerException e){ //Not nice, occurs only at startup
jbone: } // try
jbone: } else {
jbone: step = 0;
jbone: double amp = runData.getParamValue(runData.icAmplitudeNm);
jbone: double pos = runData.getParamValue(runData.icPositionNm);
jbone: double wid = runData.getParamValue(runData.icWidthNm);
jbone: double wln = runData.getParamValue(runData.icWavelengthNm);
jbone: if(icSelection.getSelectedItem().equals(BOX)){
jbone: inSolution.discretize(new ShapeBox(amp,pos,wid));
jbone: } else if(icSelection.getSelectedItem().equals(GAUSSIAN)){
jbone: inSolution.discretize(new ShapeGaussian(amp,pos,wid));
jbone: } else if(icSelection.getSelectedItem().equals(COSINUS)){
jbone: inSolution.discretize(new ShapeCosinus(amp,pos,wln));
jbone: } else if(icSelection.getSelectedItem().equals(SOLITON)){
jbone: inSolution.discretize(new ShapeSoliton(amp,pos,wid));
jbone: } else if(icSelection.getSelectedItem().equals(WAVEPACKET)){
jbone: inSolution.discretize(new ShapeWavePacket(amp,pos,wid,wln));
jbone: } else if(icSelection.getSelectedItem().equals(OPTION)){
jbone: inSolution.discretize(new ShapeOption(amp,pos,wid));
jbone: } else {
jbone: /** throw new ShapeException("Can't find shape"
jbone: + icSelection.getSelectedItem()); */
jbone: } // if
jbone: } // if
jbone: // solution.setIC(icSelection.getSelectedItem());
jbone: solution.setTime(0);
jbone: solution.rescale();
jbone: solution.previous(runData,physData);
jbone: solution.updateHeaders(runData, step);
jbone: plotArea.repaint();
jbone: } // setInitialCondition
jbone:
jbone: /** Modify defaults parameters the HTML tags from the web page
jbone: The new method introduced with V4.0 is no longuer backwards
jbone: compatible with the initial jbone V1.0 input parameters.
jbone: ****************************************************************************/
jbone: public void tagModify() {
jbone: String tmp;
jbone: if((tmp = getParameter("pde")) != null) pdeSelection.select(tmp);
jbone: if((tmp = getParameter("ic")) != null) icSelection.select(tmp);
jbone: if((tmp = getParameter("method")) != null) methodSelection.select(tmp);
jbone: if((tmp = getParameter("scheme")) != null) schemeSelection.select(tmp);
jbone: runData.tagModify(this);
jbone: } // tagModify
jbone:
jbone: /** Parameter info
jbone: ****************************************************************************/
jbone: public String[][] getParameterInfo() {
jbone: String[][] info = {
jbone: {"method", "name", "The name of the method"},
jbone: {"scheme", "name", "The name of the scheme"},
jbone: {"ic", "name", "The name of the initial condition"},
jbone: {"pde", "name", "The name of the PDE"},
jbone: {"Velocity", "double", "The advection velocity"},
jbone: {"Diffusion", "double", "The diffusion coefficient"},
jbone: {"Dispersion", "double", "The dispersion coefficient"},
jbone: {"TimeStep", "double", "The time step"},
jbone: {"MeshPoints", "int", "The number of mesh points"},
jbone: {"Walkers", "int", "The number of random walkers"},
jbone: {"TimeImplicit", "double", "Implicity parameter in time"},
jbone: {"TuneIntegr", "double", "Tunable integration parameter"},
jbone: {"ICAmplitude", "double", "Initial condition amplitude"},
jbone: {"ICPosition", "double", "Initial condition position"},
jbone: {"ICWidth", "double", "Initial condition width"},
jbone: {"Wavelength", "double", "Initial condition wavelength"},
jbone: {"RunTime", "double", "Run time"},
jbone: {"MeshLeft", "double", "Left position of the mesh"},
jbone: {"MeshLength", "double", "The length of the simulation box"},
jbone: {"PhysDataCase", "double", "Type of Physics, e.g potential barrier"},
jbone: {"PhysDataValue", "double", "Value of Physics, e.g. barrier height"}
jbone: };
jbone: return info;
jbone: } // getParameterInfo
jbone:
jbone: /** Applet start a new thread
jbone: ****************************************************************************/
jbone: public void start() {
jbone: if (runThread == null) { //Start a new thread
jbone: runThread = new Thread(this);
jbone: }
jbone: if (frozen) { //Wait for an action from the user
jbone: } else { runThread.start(); } //... or on with the calculation
jbone: }
jbone:
jbone: /** Applet stop
jbone: ****************************************************************************/
jbone: public void stop() {
jbone: runThread = null;
jbone: frozen = true;
jbone: }
jbone:
jbone: /** Contains the main loop for the time stepping.
jbone: @see ShapeFunction
jbone: @see Mesh
jbone: @see Solution
jbone: ****************************************************************************/
jbone: public void run() {
jbone: Thread.currentThread().setPriority(Thread.MIN_PRIORITY);
jbone: long startTime = System.currentTimeMillis();
jbone: Thread currentThread = Thread.currentThread();
jbone: double runTime = runData.getParamValue(runData.runTimeNm);
jbone: double timeStep = runData.getParamValue(runData.timeStepNm);
jbone:
jbone: while (currentThread == runThread && !frozen &&
jbone: ( solution.getTime() < runTime && nstep < 0 ||
jbone: nstep > 0 ) ) {
jbone: step++; nstep--; //One step forward in time
jbone: physData.update(step*timeStep);
jbone: solution.setTime(step*timeStep);
jbone: solution.updateHeaders(runData, step);
jbone: solution.next(runData, physData);
jbone:
jbone: if(operSelection.getSelectedItem().equals(CONSOLE))
jbone: solution.output(step); //Output in ASCII
jbone: plotArea.repaint(); //Update graphics
jbone:
jbone: try { startTime += delay; //Delay according to lag
jbone: Thread.sleep(Math.max(30, startTime-System.currentTimeMillis()));
jbone: } catch (InterruptedException e) { break; }
jbone: } // while
jbone: frozen=true;
jbone: } // run
jbone:
jbone: /** A new mesh is created by RunData
jbone: @see RunData
jbone: ****************************************************************************/
jbone: public void runDataNotifyMesh(){
jbone: createPhysics(); createSolution();
jbone: setInitialCondition(solution, null, false);
jbone: } // runDataNotifyMesh
jbone:
jbone: /** The number of particles is changed by RunData
jbone: @see RunData
jbone: ****************************************************************************/
jbone: public void runDataNotifyWalkers(){
jbone: createPhysics();
jbone: createSolution();
jbone: if(methodSelection.getSelectedItem().equals(PPP) ||
jbone: methodSelection.getSelectedItem().equals(PPPS) ){
jbone: Solution oldSolution = solution;
jbone: setInitialCondition(solution,new ShapeNumerical(oldSolution),true);
jbone: } else
jbone: setInitialCondition(solution,null,false);
jbone: } // runDataNotifyWalkers
jbone:
jbone: /** Responds to the user actions through the mouse and buttons (Java1.0).
jbone: Yes, we know this sucks compare to Java 1.1, but we want to be
jbone: compatible with as many browsers as possible. There is a lot of
jbone: old stuff out there...
jbone: @deprecated
jbone: ****************************************************************************/
jbone: public boolean action(Event e, Object arg) {
jbone: if(e.target instanceof Choice){
jbone: for(int i = 0; i < pdeNames.length; i++)
jbone: if(((String)arg).equals(pdeNames[i])){
jbone: solution.setPde(pdeSelection.getSelectedItem());
jbone: schemeSelection.sync(solution);
jbone: solution.setScheme(schemeSelection.getSelectedItem());
jbone: return true;
jbone: } // if
jbone: for(int i = 0; i < methodNames.length; i++)
jbone: if(((String)arg).equals(methodNames[i])){
jbone: if(methodSelection.getSelectedItem().equals(PPP)||
jbone: methodSelection.getSelectedItem().equals(PPPS) ){
jbone: Solution oldSolution = solution;
jbone: createSolution();
jbone: solution.setScheme(schemeSelection.getSelectedItem());
jbone: setInitialCondition(solution,new ShapeNumerical(oldSolution),true);
jbone: } else if (methodSelection.getSelectedItem().equals(SMP)||
jbone: methodSelection.getSelectedItem().equals(SMPS) ) {
jbone: createSolution();
jbone: solution.setScheme(schemeSelection.getSelectedItem());
jbone: setInitialCondition(solution,null,false);
jbone: } else {
jbone: Solution oldSolution = solution;
jbone: createSolution();
jbone: solution.setScheme(schemeSelection.getSelectedItem());
jbone: setInitialCondition(solution,new ShapeNumerical(oldSolution),true);
jbone: } // if
jbone: return true;
jbone: } // if
jbone: for(int i = 0; i < schemeNames.length; i++)
jbone: if(((String)arg).equals(schemeNames[i])){
jbone: Solution oldSolution = solution;
jbone: solution.setScheme(schemeSelection.getSelectedItem());
jbone: setInitialCondition(solution,new ShapeNumerical(oldSolution),true);
jbone: return true;
jbone: } // if
jbone: for(int i = 0; i < icNames.length; i++)
jbone: if(((String)arg).equals(icNames[i])){
jbone: setInitialCondition(solution, null, false);
jbone: return true;
jbone: } // if
jbone: if (EDIT.equals(operSelection.getSelectedItem()))
jbone: runData.displayTag(this);
jbone: if (EDITALL.equals(operSelection.getSelectedItem()))
jbone: runData.displayAll();
jbone: } // if
jbone: if(e.target instanceof Button) {
jbone: boolean ret = false;
jbone: if(((String)arg).equals(startName)) {
jbone: frozen = !frozen;
jbone: nstep = -1;
jbone: ret = true;
jbone: } else if(((String)arg).equals(step1Name)) {
jbone: frozen = false;
jbone: nstep = 1;
jbone: ret = true;
jbone: } else if(((String)arg).equals(step10Name)) {
jbone: frozen = false;
jbone: nstep = 10;
jbone: ret = true;
jbone: } else if(((String)arg).equals(displayName)) {
jbone: plotArea.toggleHeaders();
jbone: solution.rescale();
jbone: plotArea.repaint();
jbone: ret = true;
jbone: } else if(((String)arg).equals(initializeName)) {
jbone: frozen = true;
jbone: runData.createMesh();
jbone: createSolution();
jbone: setInitialCondition(solution, null, false);
jbone: step = 0;
jbone: nstep = -1;
jbone: ret = true;
jbone: } // if
jbone: if(ret) {
jbone: runThread = null;
jbone: start();
jbone: } // if
jbone: return ret;
jbone: } // if
jbone: return false;
jbone: } // action
jbone:
jbone: /** Destroy application
jbone: @deprecated
jbone: ****************************************************************************/
jbone: public boolean handleEvent(Event e) { //Java1.0
jbone: if (e.id == Event.WINDOW_DESTROY) { System.exit(0); }
jbone: return super.handleEvent(e);
jbone: }
jbone:
jbone: /** Print mouse coordinates to console
jbone: @deprecated
jbone: ****************************************************************************/
jbone: //Java1.0: public boolean mouseXXX(Event e, int x, int y)
jbone: //Java1.0: with XXX={Up, Down, Drag, Enter, Exit}
jbone: public boolean mouseDown(Event e, int x, int y) { //Java1.0
jbone: double[] coord = solution.measure(x,y);
jbone: double timeStep = runData.getParamValue(runData.timeStepNm);
jbone:
jbone: if (isAnApplet)
jbone: showStatus("Clicked coord ("+coord[0]+" ; "+coord[1]+")");
jbone: System.out.println("step="+step+
jbone: " time="+(step*timeStep)+
jbone: " coord=("+coord[0]+" ; "+coord[1]+")");
jbone: return true;
jbone: }
jbone:
jbone: /** Method to start the Applet as an application
jbone: @param args Not used
jbone: ****************************************************************************/
jbone: public static void main(String[] args) {
jbone: Frame f = new Frame("jbone");
jbone: jbone jBone = new jbone();
jbone: jBone.isAnApplet = false;
jbone: jBone.init();
jbone: f.add("Center", jBone);
jbone: f.pack();
jbone: f.show();
jbone: } // main
jbone:
jbone: } // jbone
BandMatrix: /* $Id: BandMatrix.java,v 1.9 2000/05/26 11:48:23 andrej Exp $ */
BandMatrix:
BandMatrix: /***************************************************************************
BandMatrix: * BandMatrix - Linear algebra package for band matrices
BandMatrix: ***************************************************************************/
BandMatrix: public class BandMatrix {
BandMatrix: double[][] m; //Matrix
BandMatrix: double[] d; //Temporary storage for diagonal
BandMatrix: double det;
BandMatrix: int diagos, lines, n;
BandMatrix: boolean isDecomposed = false;
BandMatrix: int l=0; //Band indices (left, center, right)
BandMatrix: int c=1;
BandMatrix: int r=2;
BandMatrix:
BandMatrix: /** Constructor */
BandMatrix: public BandMatrix(int diagos, int lines) {
BandMatrix: m = new double[diagos][lines];
BandMatrix: this.diagos =diagos;
BandMatrix: this.lines =lines;
BandMatrix: this.n =lines-1;
BandMatrix: }
BandMatrix:
BandMatrix: /** Store matrix elements */
BandMatrix: public void set(int i, int j, double v) { m[j][i]=v; }
BandMatrix: public void setL(int i, double v) { m[l][i]=v; }
BandMatrix: public void setD(int i, double v) { m[c][i]=v; }
BandMatrix: public void setR(int i, double v) { m[r][i]=v; }
BandMatrix:
BandMatrix: /** Retrieve value of matrix elements */
BandMatrix: public double get(int j, int i) { return m[j][i]; }
BandMatrix: public double getL(int i) { return m[l][i]; }
BandMatrix: public double getD(int i) { return m[c][i]; }
BandMatrix: public double getR(int i) { return m[r][i]; }
BandMatrix:
BandMatrix: /** Matrix times vector */
BandMatrix: public double[] dot(double[] v) {
BandMatrix: double[] res = new double[lines];
BandMatrix: int hw=(diagos-1)/2;
BandMatrix: for (int i=0; i<=n; i++) {
BandMatrix: for (int j=Math.max(0,hw-i); j<=Math.min(2,hw+n-i); j++) {
BandMatrix: res[i]=res[i]+m[j][i]*v[i+j-1];
BandMatrix: }
BandMatrix: }
BandMatrix: return res;
BandMatrix: }
BandMatrix:
BandMatrix: /** Direct LU solver for 3-banded matrix with multiple RHS and BC
BandMatrix: The matrix is decomposed once only for the first rhs vector
BandMatrix: Periodic boundary conditions are taken care of by the upper-left
BandMatrix: lower-right most matrix elements
BandMatrix: @param rhs Right hand side vector
BandMatrix: @return Solution of the linear system
BandMatrix: ****************************************************************************/
BandMatrix: public double[] solve3(double[] rhs) {
BandMatrix: double[] sol = new double[lines];
BandMatrix: int i;
BandMatrix:
BandMatrix: if (!isDecomposed) { //Forward substitution of matrix
BandMatrix: d = new double[lines];
BandMatrix: d[0]= m[c][0];
BandMatrix: m[r][0]= m[r][0]/d[0];
BandMatrix: m[l][0]=-m[l][0]/d[0];
BandMatrix: for (i=1; i<n; i++) {
BandMatrix: d[i] = m[c][i]-m[l][i]*m[r][i-1];
BandMatrix: m[r][i]= m[r][i]/d[i];
BandMatrix: m[c][i]=-m[l][i]/d[i];
BandMatrix: m[l][i]=-m[l][i]/d[i]*m[l][i-1];
BandMatrix: }
BandMatrix: d[n]=m[c][n]+m[l][n]*(m[l][n-1]-m[r][n-1]);
BandMatrix: m[r][n]= m[r][n]/d[n];
BandMatrix: m[c][n]=-m[l][n]/d[n];
BandMatrix:
BandMatrix: m[l][n-1]=m[l][n-1]-m[r][n-1]; //Backward substitution of matrix
BandMatrix: for (i=n-2; i>=0; i--) {
BandMatrix: m[l][i]= m[l][i]-m[r][i]*m[l][i+1];
BandMatrix: }
BandMatrix: det = 1.0 + m[r][n]*m[l][0];
BandMatrix: isDecomposed=true;
BandMatrix: }
BandMatrix:
BandMatrix: rhs[0]=rhs[0]/d[0];
BandMatrix: for (i=1; i<=n; i++) { //Forward substitution of rhs
BandMatrix: rhs[i]=rhs[i]/d[i] +m[c][i]*rhs[i-1];
BandMatrix: }
BandMatrix: for (i=n-2; i>=0; i--) { //Backward substitution of rhs
BandMatrix: rhs[i]=rhs[i] -m[r][i]*rhs[i+1];
BandMatrix: }
BandMatrix:
BandMatrix: sol[0]=(rhs[0]+m[l][0]*rhs[n])/det; //Built-in periodicity
BandMatrix: sol[n]=(rhs[n]-m[r][n]*rhs[0])/det;
BandMatrix: for (i=1;i<n;i++) {
BandMatrix: sol[i]=rhs[i]+m[l][i]*sol[n];
BandMatrix: }
BandMatrix: return sol;
BandMatrix: }
BandMatrix:
BandMatrix:
BandMatrix: /** Projected Symmetric Over Relaxation for 3-banded matrix
BandMatrix: Periodic boundary conditions are taken care of by the upper-left
BandMatrix: lower-right most matrix elements
BandMatrix: @param rhs Right hand side vector
BandMatrix: @param ini Initial iterate
BandMatrix: @param min Lower limit for obstacle problems
BandMatrix: @param max Upper limit for obstacle problems
BandMatrix: @param precision Relative precision required
BandMatrix: @param w Relaxation parameter to be chosen within [0;2]
BandMatrix: @param maxIterations Maximum number of iterations
BandMatrix: @return Solution of the linear system
BandMatrix: ****************************************************************************/
BandMatrix: public double[] ssor3(double[] rhs,double[] ini,double[] min,double[] max,
BandMatrix: double precision,double w,int maxIterations) {
BandMatrix: double[] hlf = new double[lines]; //Half-step
BandMatrix: double[] sol = new double[lines]; //Solution
BandMatrix: boolean converged = false;
BandMatrix: int iterations = 0;
BandMatrix: int i;
BandMatrix: while (!converged && iterations<maxIterations){
BandMatrix:
BandMatrix: //Gauss Seidel, forward one half-step
BandMatrix: hlf[0]=Math.max(min[0],Math.min(max[0],
BandMatrix: ( w *(rhs[0]-m[r][0]*ini[1]-m[l][0]*hlf[n])+
BandMatrix: (1-w)*m[c][0]*ini[0])/m[c][0] ));
BandMatrix: for(i=1;i<=n-1;i++)
BandMatrix: hlf[i]=Math.max(min[i],Math.min(max[i],
BandMatrix: ( w *(rhs[i]-m[r][i]*ini[i+1]-m[l][i]*hlf[i-1])+
BandMatrix: (1-w)*m[c][i]*ini[i])/m[c][i] ));
BandMatrix: hlf[n]=Math.max(min[n],Math.min(max[n],
BandMatrix: ( w *(rhs[n]-m[r][n]*ini[0]-m[l][n]*hlf[n-1])+
BandMatrix: (1-w)*m[c][n]*ini[n])/m[c][n] ));
BandMatrix:
BandMatrix: //Gauss Seidel, backward one half-step
BandMatrix: sol[n]=Math.max(min[n],Math.min(max[n],
BandMatrix: ( w *(rhs[n]-m[l][n]*hlf[n-1]-m[r][n]*sol[0])+
BandMatrix: (1-w)*m[c][n]*hlf[n])/m[c][n] ));
BandMatrix: for(i=n-1;i>=1;i--)
BandMatrix: sol[i]=Math.max(min[i],Math.min(max[i],
BandMatrix: ( w *(rhs[i]-m[l][i]*hlf[i-1]-m[r][i]*sol[i+1])+
BandMatrix: (1-w)*m[c][i]*hlf[i])/m[c][i] ));
BandMatrix: sol[0]=Math.max(min[0],Math.min(max[0],
BandMatrix: ( w *(rhs[0]-m[l][0]*hlf[n]-m[r][0]*sol[1])+
BandMatrix: (1-w)*m[c][0]*hlf[0])/m[c][0] ));
BandMatrix:
BandMatrix: converged=isEqual(sol,ini,precision);
BandMatrix: iterations++;
BandMatrix:
BandMatrix: //Diagnostic
BandMatrix: if (false && (iterations<30 || converged) ) {
BandMatrix: double err=0.;
BandMatrix: for (i=0;i<=n;i++) {
BandMatrix: err=Math.max(err,Math.abs(sol[i]-ini[i]));
BandMatrix: }
BandMatrix: System.out.println("PSSOR iteration "+iterations+", max(err)="+err);
BandMatrix: }
BandMatrix: for (i=0;i<=n;i++) {ini[i]=sol[i];}
BandMatrix: }
BandMatrix: return sol;
BandMatrix:
BandMatrix: }//ssor3
BandMatrix:
BandMatrix:
BandMatrix: /** Projected Symmetric Over Relaxation for 3-banded matrix
BandMatrix: Periodic boundary conditions are taken care of by the upper-left
BandMatrix: lower-right most matrix elements
BandMatrix: @param rhs Right hand side vector
BandMatrix: @param ini Initial iterate
BandMatrix: @return Solution of the linear system
BandMatrix: ****************************************************************************/
BandMatrix: public double[] ssor3(double[] rhs,double[] ini) {
BandMatrix: double precision = Math.pow(10.,-5); //Default SSOR parameters
BandMatrix: int maxIterations= 1000;
BandMatrix: double w = 1.2;
BandMatrix: double[] min = new double[n+1]; //No projection limits
BandMatrix: double[] max = new double[n+1];
BandMatrix: for (int i=0; i<=n; i++) {
BandMatrix: min[i]=Double.NEGATIVE_INFINITY;
BandMatrix: max[i]=Double.POSITIVE_INFINITY;
BandMatrix: }
BandMatrix: return ssor3(rhs, ini, min, max, precision, w, maxIterations);
BandMatrix: }
BandMatrix:
BandMatrix:
BandMatrix:
BandMatrix: /** Decide whether two vectors are equal to a certain degree of precision
BandMatrix: @param sol First vector
BandMatrix: @param ini Second vector
BandMatrix: @param precision Absolute precision
BandMatrix: @return true is both vectors are equal
BandMatrix: ****************************************************************************/
BandMatrix: public boolean isEqual(double[] sol, double[] ini, double precision) {
BandMatrix: test: for (int k=0 ; k<=n ; k++){
BandMatrix: if (Math.abs(sol[k]-ini[k]) >= precision )
BandMatrix: { return false; } else { continue test; }
BandMatrix: }
BandMatrix: return true;
BandMatrix: }//isEqual
BandMatrix:
BandMatrix: } //End of BandMatrix
BandMatrixC: /* $Id: BandMatrixC.java,v 1.2 2000/05/05 17:32:34 andrej Exp $ */
BandMatrixC:
BandMatrixC: /***************************************************************************
BandMatrixC: * BandMatrixC - Linear algebra package for Complex band matrices
BandMatrixC: ***************************************************************************/
BandMatrixC: public class BandMatrixC {
BandMatrixC: Complex[][] m; //Matrix
BandMatrixC: Complex[] d; //Temporary storage for diagonal
BandMatrixC: Complex det;
BandMatrixC: int diagos, lines, n;
BandMatrixC: boolean isDecomposed = false;
BandMatrixC: int l=0; //Band indices (left, center, right)
BandMatrixC: int c=1;
BandMatrixC: int r=2;
BandMatrixC:
BandMatrixC: /** Constructor */
BandMatrixC: public BandMatrixC(int diagos, int lines) {
BandMatrixC: m = new Complex[diagos][lines];
BandMatrixC: this.diagos =diagos;
BandMatrixC: this.lines =lines;
BandMatrixC: this.n =lines-1;
BandMatrixC: }
BandMatrixC:
BandMatrixC: /** Store matrix elements */
BandMatrixC: public void set(int i, int j, Complex v) { m[j][i]=v; }
BandMatrixC: public void setL(int i, Complex v) { m[l][i]=v; }
BandMatrixC: public void setD(int i, Complex v) { m[c][i]=v; }
BandMatrixC: public void setR(int i, Complex v) { m[r][i]=v; }
BandMatrixC:
BandMatrixC: /** Retrieve value of matrix elements */
BandMatrixC: public Complex get(int j, int i) { return m[j][i]; }
BandMatrixC: public Complex getL(int i) { return m[l][i]; }
BandMatrixC: public Complex getD(int i) { return m[c][i]; }
BandMatrixC: public Complex getR(int i) { return m[r][i]; }
BandMatrixC:
BandMatrixC: /** Matrix times vector */
BandMatrixC: public Complex[] dot(Complex[] v) {
BandMatrixC: Complex[] res = new Complex[this.lines];
BandMatrixC: Complex z = new Complex();
BandMatrixC: int hw=(diagos-1)/2;
BandMatrixC: for (int i=0; i<=n; i++) {
BandMatrixC: res[i]=new Complex(0.);
BandMatrixC: for (int j=Math.max(0,hw-i); j<=Math.min(2,hw+n-i); j++) {
BandMatrixC: z =m[j][i].mul(v[i+j-1]);
BandMatrixC: res[i]=res[i].add(z);
BandMatrixC: }
BandMatrixC: }
BandMatrixC: return res;
BandMatrixC: }
BandMatrixC:
BandMatrixC: /** Gaussian elimination for 3-banded matrix with multiple RHS and BC */
BandMatrixC: public Complex[] solve3(Complex[] rhs) {
BandMatrixC: Complex[] sol = new Complex[lines];
BandMatrixC: Complex z1 = new Complex();
BandMatrixC: Complex z2 = new Complex();
BandMatrixC: int i;
BandMatrixC:
BandMatrixC: if (!isDecomposed) { //Forward substitution of matrix
BandMatrixC: d = new Complex[lines];
BandMatrixC: d[0]= new Complex(m[c][0]);
BandMatrixC: m[r][0]= m[r][0].div(d[0]);
BandMatrixC: m[l][0]=(new Complex(-1.0)).mul(m[l][0].div(d[0]));
BandMatrixC:
BandMatrixC: for (i=1; i<n; i++) {
BandMatrixC: d[i]=new Complex();
BandMatrixC: z1 = m[l][i].mul(m[r][i-1]);
BandMatrixC: d[i]=m[c][i].sub(z1);
BandMatrixC: m[r][i]= m[r][i].div(d[i]);
BandMatrixC: m[c][i]=(new Complex(-1.0)).mul(m[l][i].div(d[i]));
BandMatrixC: z1=(new Complex(-1.0)).mul(m[l][i].div(d[i]));
BandMatrixC: m[l][i]=z1.mul(m[l][i-1]);
BandMatrixC: }
BandMatrixC: z1=m[l][n-1].sub(m[r][n-1]);
BandMatrixC: z2=m[l][n].mul(z1);
BandMatrixC: d[n]=new Complex();
BandMatrixC: d[n]=m[c][n].add(z2);
BandMatrixC: m[r][n]= m[r][n].div(d[n]);
BandMatrixC: m[c][n]=(new Complex(-1.0)).mul(m[l][n].div(d[n]));
BandMatrixC:
BandMatrixC: m[l][n-1]=m[l][n-1].sub(m[r][n-1]); //Backward substitution of matrix
BandMatrixC:
BandMatrixC: for (i=n-2; i>=0; i--) {
BandMatrixC: z1=m[r][i].mul(m[l][i+1]);
BandMatrixC: m[l][i]=m[l][i].sub(z1);
BandMatrixC: }
BandMatrixC: z1=m[r][n].mul(m[l][0]);
BandMatrixC: det=(new Complex(1.0)).add(z1);
BandMatrixC: isDecomposed=true;
BandMatrixC: }
BandMatrixC: rhs[0]=rhs[0].div(d[0]);
BandMatrixC:
BandMatrixC: for (i=1; i<=n; i++) { //Forward substitution of rhs
BandMatrixC: z1=m[c][i].mul(rhs[i-1]);
BandMatrixC: z2=rhs[i].div(d[i]);
BandMatrixC: rhs[i]=z2.add(z1);
BandMatrixC: }
BandMatrixC: for (i=n-2; i>=0; i--) { //Backward substitution of rhs
BandMatrixC: z1=m[r][i].mul(rhs[i+1]);
BandMatrixC: rhs[i]=rhs[i].sub(z1);
BandMatrixC: }
BandMatrixC: z1=m[l][0].mul(rhs[n]);
BandMatrixC: z2=rhs[0].add(z1);
BandMatrixC: sol[0]=z2.div(det);
BandMatrixC: //Built-in periodicity
BandMatrixC: z1=m[r][n].mul(rhs[0]);
BandMatrixC: z2=rhs[n].sub(z1);
BandMatrixC: sol[n]=z2.div(det);
BandMatrixC: for (i=1;i<n;i++) {
BandMatrixC: z1=m[l][i].mul(sol[n]);
BandMatrixC: sol[i]=rhs[i].add(z1);
BandMatrixC: }
BandMatrixC: return sol;
BandMatrixC: }
BandMatrixC:
BandMatrixC: } //End of BandMatrixC
CHASolution: /* $Id: CHASolution.java,v 1.26 2001/01/15 14:22:21 andrej Exp $ */
CHASolution:
CHASolution: /** The class CHASolution contains the Characteristics solver.
CHASolution: @version $Revision: 1.26 $
CHASolution: */
CHASolution: public class CHASolution extends FluidSolution{
CHASolution: /** Creates a CHASolution object.
CHASolution: @param runData The run time parameters
CHASolution: */
CHASolution: public CHASolution(RunData runData){
CHASolution: super(runData);
CHASolution: } // FDSolution
CHASolution:
CHASolution:
CHASolution: /** Advance the solution forward one step in time.
CHASolution: The equation is set by the internal variable pde and
CHASolution: the numerical scheme by the internal variable scheme
CHASolution: @param runData List of run parameters
CHASolution: @see RunData
CHASolution: @return True if a scheme is implemented for that equation
CHASolution: ****************************************************************************/
CHASolution: public boolean next(RunData runData, PhysData physData) {
CHASolution: int n=f.length-1;
CHASolution: boolean isDefined=true;
CHASolution:
CHASolution: for (int j=0; j<=n; j++) { gp[j]=g[j]; } //Plot initial condition
CHASolution:
CHASolution: if(pde.equals(jbone.ADVECTION) &&
CHASolution: scheme.equals(jbone.CUBFEM)){
CHASolution: //=======================================================================
CHASolution: // CHA - Advection/diffusion - CubicFEM - Yabe+Aoki CPC 66(1991)219
CHASolution: //=======================================================================
CHASolution: double timeStep= runData.getParamValue(runData.timeStepNm);
CHASolution: double velocity= runData.getParamValue(runData.velocityNm);
CHASolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
CHASolution:
CHASolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
CHASolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
CHASolution: double a,b;
CHASolution:
CHASolution: for (int j=0; j<n; j++) {
CHASolution: a=dx[0]*(df[j]+ df[j+1])-2*(f[j+1]-f[j]);
CHASolution: b=dx[0]*(df[j]+2*df[j+1])-3*(f[j+1]-f[j]);
CHASolution: fp[j+1]= f[j+1] -beta*(dx[0]*df[j+1]-beta*(b-beta*a));
CHASolution: dfp[j+1]= df[j+1] -beta/dx[0]*(2*b-3*beta*a);
CHASolution: }
CHASolution: a=dx[0]*(df[n]+ df[0])-2*(f[0]-f[n]);
CHASolution: b=dx[0]*(df[n]+2*df[0])-3*(f[0]-f[n]);
CHASolution: fp[0]= f[0] -beta*(dx[0]*df[0]-beta*(b-beta*a));
CHASolution: dfp[0]= df[0] -beta/dx[0]*(2*b-3*beta*a);
CHASolution:
CHASolution: } else if(pde.equals(jbone.ADVECTION) &&
CHASolution: scheme.equals(jbone.CUBSPL)){
CHASolution: //=======================================================================
CHASolution: // CHA - Advection/diffusion - Cubic splines
CHASolution: //=======================================================================
CHASolution: double timeStep= runData.getParamValue(runData.timeStepNm);
CHASolution: double velocity= runData.getParamValue(runData.velocityNm);
CHASolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
CHASolution:
CHASolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
CHASolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
CHASolution: double h = dx[0];
CHASolution: BandMatrix a= new BandMatrix(3, f.length);
CHASolution: double[] c= new double[f.length];
CHASolution: int j, j1, j2;
CHASolution: double z, z1, z2, dm, dp;
CHASolution:
CHASolution: for (j=0; j<=n; j++) { // Matrix
CHASolution: a.setL(j, h/6 );
CHASolution: a.setD(j, 4*h/6 );
CHASolution: a.setR(j, h/6 );
CHASolution: }
CHASolution:
CHASolution: dp= (f[0]-f[n])/h; //Right hand side 2nd derivative
CHASolution: for (j=0; j<n; j++) {
CHASolution: dm= dp; dp= (f[j+1]-f[j])/h;
CHASolution: c[j]= dp-dm;
CHASolution: }
CHASolution: dm= dp; dp= (f[0]-f[n])/h;
CHASolution: c[n]= dp-dm;
CHASolution:
CHASolution: dfp=a.solve3(c); //Solve linear problem
CHASolution:
CHASolution: double boxLen= (double)mesh.size() * dx[0];
CHASolution:
CHASolution: for (j=0; j<=n; j++) { //Propagate along characteristics
CHASolution: z= x[j]-timeStep*velocity; //Allow arbitrary time step
CHASolution: while (z> boxLen) { z= z-boxLen* (int)( z/boxLen); }
CHASolution: while (z< 0) { z= z+boxLen*(1+(int)(-z/boxLen)); }
CHASolution: j1= (int)(z/h); j2=j1+1;
CHASolution: if (j1 == n) { j2= 0; }
CHASolution: z1 = (z-x[j1])/h; z2 = 1-z1; //Cubic spline interpolation
CHASolution: fp[j] = z2*f[j1] +z1*f[j2] +(h*h)/6*((z2*z2*z2 -z2)*dfp[j1] +
CHASolution: (z1*z1*z1 -z1)*dfp[j2] );
CHASolution: }
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.EXC6_01)){
CHASolution: //=======================================================================
CHASolution: // CHA - Exercise 6.01
CHASolution: //=======================================================================
CHASolution: //INCLUDE_sol6.01.java_INCLUDE//
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.EXC6_02)){
CHASolution: //=======================================================================
CHASolution: // CHA - Exercise 6.02
CHASolution: //=======================================================================
CHASolution: //INCLUDE_sol6.02.java_INCLUDE//
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.EXC6_03)){
CHASolution: //=======================================================================
CHASolution: // CHA - Exercise 6.03
CHASolution: //=======================================================================
CHASolution: //INCLUDE_sol6.03.java_INCLUDE//
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.EXC6_04)){
CHASolution: //=======================================================================
CHASolution: // CHA - Exercise 6.04
CHASolution: //=======================================================================
CHASolution: //INCLUDE_sol6.04.java_INCLUDE//
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.EXC6_05)){
CHASolution: //=======================================================================
CHASolution: // CHA - Exercise 6.05
CHASolution: //=======================================================================
CHASolution: //INCLUDE_sol6.05.java_INCLUDE//
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.EXCF_01)){
CHASolution: //=======================================================================
CHASolution: // CHA - Exercise F.01
CHASolution: //=======================================================================
CHASolution: //INCLUDE_solF.01.java_INCLUDE//
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.MY_SCHEME)){
CHASolution: //=======================================================================
CHASolution: // CHA - Project 1
CHASolution: //=======================================================================
CHASolution: //INCLUDE_prj6.01.java_INCLUDE//
CHASolution:
CHASolution: } else if(pde.equals(jbone.EXERCISE)
CHASolution: && scheme.equals(jbone.MY_SCHEME_2)){
CHASolution: //=======================================================================
CHASolution: // CHA - Project 2
CHASolution: //=======================================================================
CHASolution: //INCLUDE_prj6.02.java_INCLUDE//
CHASolution:
CHASolution: } // if
CHASolution: if(isDefined)
CHASolution: shiftLevels();
CHASolution: return isDefined;
CHASolution: } // next
CHASolution:
CHASolution: /** Take the solution backward one step to initialize schemes
CHASolution: with 3 time levels.
CHASolution: The equation is set by the internal variable pde and
CHASolution: the numerical scheme by the internal variable scheme
CHASolution: @param runData List of run parameters
CHASolution: @see RunData
CHASolution: @return True if an initialisation scheme is implemented for that equation
CHASolution: ****************************************************************************/
CHASolution: public boolean previous(RunData runData, PhysData physData){
CHASolution: int n=f.length-1;
CHASolution: boolean isDefined=false;
CHASolution:
CHASolution: if(pde.equals(jbone.ADVECTION) &&
CHASolution: scheme.equals(jbone.CUBFEM)){
CHASolution: //=======================================================================
CHASolution: // CHA - Advection/diffusion - CubicFEM - Yabe+Aoki CPC 66(1991)219
CHASolution: //=======================================================================
CHASolution: for (int i=1; i<n; i++) {
CHASolution: df[i] =(f[i+1]-f[i-1])/(2.*dx[0]);
CHASolution: }
CHASolution: df[0]=(f[1]-f[n ])/(2.*dx[0]);
CHASolution: df[n]=(f[0]-f[n-1])/(2.*dx[0]);
CHASolution: isDefined=true;
CHASolution: }
CHASolution: return isDefined;
CHASolution: } // previous
CHASolution:
CHASolution: /** Tells whether the solution implements a option
CHASolution: @param option The the option to implement
CHASolution: @return True if the option is implemented
CHASolution: @see jbone#pdeNames
CHASolution: @see jbone#schemeNames
CHASolution: */
CHASolution: public boolean hasOption(String option){
CHASolution: try {
CHASolution: if(option.equals(jbone.ADVECTION))
CHASolution: return true;
CHASolution: else if (option.equals(jbone.EXERCISE))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.CUBFEM))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.CUBSPL))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC6_01))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC6_02))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC6_03))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC6_04))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC6_05))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXCF_01))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME))
CHASolution: return true;
CHASolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME_2))
CHASolution: return true;
CHASolution: return false;
CHASolution: } catch (NullPointerException e){
CHASolution: return false;
CHASolution: } // try
CHASolution: } // hasOption
CHASolution:
CHASolution: } // class CHASolution
FDSolution: /* $Id: FDSolution.java,v 1.35 2001/01/16 18:35:55 andrej Exp $ */
FDSolution:
FDSolution: /******************************************************************************
FDSolution: * FDSolution -- contains the finite difference solver.
FDSolution: * @see FluidSolution
FDSolution: * @see Solution
FDSolution: ******************************************************************************/
FDSolution: public class FDSolution extends FluidSolution{
FDSolution:
FDSolution: /** Creates a FDSolution object.
FDSolution: @param runData The run time parameters
FDSolution: */
FDSolution: public FDSolution(RunData runData){
FDSolution: super(runData);
FDSolution: } // FDSolution
FDSolution:
FDSolution: /** Advance the solution forward one step in time.
FDSolution: The equation is set by the internal variable pde and
FDSolution: the numerical scheme by the internal variable scheme
FDSolution: @param runData List of run parameters
FDSolution: @see RunData
FDSolution: @return True if a scheme is implemented for that equation
FDSolution: ****************************************************************************/
FDSolution: public boolean next(RunData runData, PhysData physData) {
FDSolution: int n=f.length-1;
FDSolution: boolean isDefined=true;
FDSolution:
FDSolution: for (int j=0; j<=n; j++) { gp[j]=g[j]; } //Default initial condition
FDSolution:
FDSolution: if(pde.equals(jbone.ADVECTION)
FDSolution: && scheme.equals(jbone.EXP2LVL)){
FDSolution: //=======================================================================
FDSolution: // FD - Advection/diffusion - explicit 2 levels (Godunov advection)
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FDSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FDSolution:
FDSolution: for (int j=1; j<n; j++) {
FDSolution: fp[j]=f[j] -beta *(f[j]-f[j-1])+alpha*(f[j+1]-2.*f[j]+f[j-1]); }
FDSolution: fp[0]=f[0] -beta *(f[0]-f[ n ])+alpha*(f[ 1 ]-2.*f[0]+f[ n ]);
FDSolution: fp[n]=f[n] -beta *(f[n]-f[n-1])+alpha*(f[ 0 ]-2.*f[n]+f[n-1]);
FDSolution:
FDSolution: } else if(pde.equals(jbone.ADVECTION)
FDSolution: && scheme.equals(jbone.EXP3LVL)){
FDSolution: //=======================================================================
FDSolution: // FD - Advection/diffusion - explicit 3 levels
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FDSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FDSolution:
FDSolution: for (int j=1; j<n; j++) {
FDSolution: fp[j]=fm[j] -beta*(f[j+1]-f[j-1]) +2*alpha*(f[j+1]-2.*f[j]+f[j-1]);}
FDSolution: fp[0]=fm[0] -beta*(f[ 1 ]-f[ n ]) +2*alpha*(f[ 1 ]-2.*f[0]+f[ n ]);
FDSolution: fp[n]=fm[n] -beta*(f[ 0 ]-f[n-1]) +2*alpha*(f[ 0 ]-2.*f[n]+f[n-1]);
FDSolution:
FDSolution: } else if(pde.equals(jbone.ADVECTION)
FDSolution: && scheme.equals(jbone.IMP2LVL)){
FDSolution: //=======================================================================
FDSolution: // FD - Advection/diffusion - implicit (Crank-Nicholson diffusion)
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FDSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FDSolution: BandMatrix a = new BandMatrix(3, f.length);
FDSolution: BandMatrix b = new BandMatrix(3, f.length);
FDSolution: double[] c = new double[f.length];
FDSolution:
FDSolution: for (int j=0; j<=n; j++) {
FDSolution: a.setL(j,-0.25*beta -0.5*alpha); //Matrix elements
FDSolution: a.setD(j, 1. +alpha);
FDSolution: a.setR(j, 0.25*beta -0.5*alpha);
FDSolution: b.setL(j, 0.25*beta +0.5*alpha); //Right hand side
FDSolution: b.setD(j, 1. -alpha);
FDSolution: b.setR(j,-0.25*beta +0.5*alpha);
FDSolution: }
FDSolution:
FDSolution: c=b.dot(f); //Right hand side
FDSolution: c[0]=c[0]+b.getL(0)*f[n]; // with periodicity
FDSolution: c[n]=c[n]+b.getR(n)*f[0];
FDSolution:
FDSolution: fp=a.solve3(c); //Solve linear problem
FDSolution:
FDSolution: } else if(pde.equals(jbone.ADVECTION)
FDSolution: && scheme.equals(jbone.EXPLAX)){
FDSolution: //=======================================================================
FDSolution: // FD - Advection - explicit Lax-Wendroff - precision in O(delta t^3)
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FDSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FDSolution:
FDSolution: for (int j=1; j<n; j++) {
FDSolution: fp[j]=f[j] -0.5*beta *(f[j+1]-f[j-1])
FDSolution: +0.5*beta*beta*(f[j+1]-2.*f[j]+f[j-1]); }
FDSolution: fp[0]=f[0] -0.5*beta *(f[ 1 ]-f[ n ])
FDSolution: +0.5*beta*beta*(f[ 1 ]-2.*f[0]+f[ n ]);
FDSolution: fp[n]=f[n] -0.5*beta *(f[ 0 ]-f[n-1])
FDSolution: +0.5*beta*beta*(f[ 0 ]-2.*f[n]+f[n-1]);
FDSolution:
FDSolution: } else if(pde.equals(jbone.ADVECTION)
FDSolution: && scheme.equals(jbone.IMPLAX)){
FDSolution: //=======================================================================
FDSolution: // FD - Advection - implicit Lax-Wendroff
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FDSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FDSolution:
FDSolution: BandMatrix a = new BandMatrix(3, f.length);
FDSolution: BandMatrix b = new BandMatrix(3, f.length);
FDSolution: double[] c = new double[f.length];
FDSolution:
FDSolution: for (int j=0; j<=n; j++) {
FDSolution: a.setL(j, (1.-beta)/(1.+beta)); //Matrix elements
FDSolution: a.setD(j, 1.);
FDSolution: a.setR(j, 0.);
FDSolution: b.setL(j, 1.); //Right hand side
FDSolution: b.setD(j, (1.-beta)/(1.+beta));
FDSolution: b.setR(j, 0.);
FDSolution: }
FDSolution:
FDSolution: c=b.dot(f); //Right hand side
FDSolution: c[0]=c[0]+b.getL(0)*f[n]; // with periodicity
FDSolution: c[n]=c[n]+b.getR(n)*f[0];
FDSolution:
FDSolution: fp=a.solve3(c); //Solve linear problem
FDSolution:
FDSolution: } else if(pde.equals(jbone.ADVECTION)
FDSolution: && scheme.equals(jbone.LAXFRIED)){
FDSolution: //=======================================================================
FDSolution: // FD - Advection - Lax-Friedrichs
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FDSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FDSolution:
FDSolution: for (int i=1; i<n; i++) {
FDSolution: fp[i]=0.5*(f[i+1]+f[i-1])
FDSolution: -0.5*beta*(f[i+1]-f[i-1]) +2*alpha*(f[i+1]-2.*f[i]+f[i-1]);}
FDSolution: fp[0]=0.5*(f[1]+f[n ])-0.5*beta*(f[1]-f[n]) +2*alpha*(f[1]-2.*f[0]+f[n]);
FDSolution: fp[n]=0.5*(f[0]+f[n-1])-0.5*beta*(f[0]-f[n-1])+2*alpha*(f[0]-2.*f[n]+f[n-1]);
FDSolution:
FDSolution: } else if(pde.equals(jbone.ADVECTION)
FDSolution: && scheme.equals(jbone.LPFROG)){
FDSolution: //=======================================================================
FDSolution: // FDTD - Advection/Wave - Leap-frog
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FDSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FDSolution:
FDSolution: for (int j=1; j<=n; j++) { //time + time_step
FDSolution: fp[j]=f[j] -beta*(g[j]-g[j-1]); }
FDSolution: fp[0]=f[0] -beta*(g[0]-g[n]);
FDSolution:
FDSolution: for (int j=0; j<=n-1; j++) { //time + 1.5*time_step
FDSolution: gp[j]=g[j] -beta*(fp[j+1]-fp[j]); }
FDSolution: gp[n]=g[n] -beta*(fp[0]-fp[n]);
FDSolution:
FDSolution: } else if(pde.equals(jbone.WAVE)
FDSolution: && scheme.equals(jbone.FDTDABS)){
FDSolution: //=======================================================================
FDSolution: // FDTD - Wave - Refraction and absobing boundary conditions
FDSolution: //=======================================================================
FDSolution:
FDSolution: } else if(pde.equals(jbone.BURGER)
FDSolution: && scheme.equals(jbone.EXP2LVL)){
FDSolution: //=======================================================================
FDSolution: // FD - Burger's shock wave equation - explicit 2 levels
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution:
FDSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);// Homogeneous problem
FDSolution:
FDSolution: double fm1 = f[ n ]; //Extra mesh points for periodicity
FDSolution: double f0 = f[ 0 ];
FDSolution: double fp1 = f[ 1 ];
FDSolution:
FDSolution: for (int j=0; j<=n; ++j) {
FDSolution: fp[j]= f[j]
FDSolution: -0.5*timeStep/dx[0]* f0*(fp1-fm1) //Wave breaking term
FDSolution: +alpha*(fp1-2*f0+fm1); //Diffusive term
FDSolution: fm1=f0; f0 =fp1;
FDSolution: fp1=f[(j+2) % (n+1)];
FDSolution: }
FDSolution:
FDSolution: } else if(pde.equals(jbone.KDV)
FDSolution: && scheme.equals(jbone.EXP3LVL)){
FDSolution: //=======================================================================
FDSolution: // FD - KdV - Scheme Zabusky & Kruskal, Phys.Rev.Lett. 15(1965)240
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double disperCo= runData.getParamValue(runData.dispersionNm);
FDSolution:
FDSolution: double nolin=timeStep/(dx[0]);
FDSolution: double gamma=timeStep/(dx[0]*dx[0]*dx[0])*disperCo;
FDSolution:
FDSolution: double fm2 = f[n-1]; double fm1 = f[ n ];
FDSolution: double f0 = f[ 0 ]; double fp1 = f[ 1 ]; double fp2 = f[ 2 ];
FDSolution:
FDSolution: for (int j=0; j<=n; ++j) {
FDSolution: fp[j]= fm[j]
FDSolution: -nolin*(fp1+f0+fm1)/3.*(fp1-fm1) //Wave breaking term
FDSolution: -gamma*(fp2 -2.*(fp1-fm1) -fm2); //Dispersive term
FDSolution:
FDSolution: fm2=fm1; fm1=f0; f0 =fp1; fp1=fp2; //Periodicity
FDSolution: fp2=f[(j+3) % (n+1)];
FDSolution: }
FDSolution:
FDSolution: } else if(pde.equals(jbone.SCHROED)
FDSolution: && scheme.equals(jbone.IMP2LVL)){
FDSolution: //=======================================================================
FDSolution: // FD - Schroedinger - implicit (Crank-Nicholson) - Cayley's form
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution:
FDSolution: BandMatrixC a = new BandMatrixC(3, h.length);
FDSolution: BandMatrixC b = new BandMatrixC(3, h.length);
FDSolution: Complex[] c = new Complex[h.length];
FDSolution: Complex z = new Complex();
FDSolution:
FDSolution: double[] V = physData.getPotential(runData); //Here Heavyside(x-L/2)
FDSolution: double scale = 10; // times scaling factor
FDSolution: double dtodx2 = timeStep/(dx[0]*dx[0]);
FDSolution: Complex pjh = new Complex(0., 0.5*dtodx2);
FDSolution: Complex mjh = new Complex(0.,-0.5*dtodx2);
FDSolution: Complex pjp1 = new Complex(1., dtodx2);
FDSolution: Complex mjp1 = new Complex(1., -dtodx2);
FDSolution:
FDSolution: for (int j=0; j<=n; j++) {
FDSolution: z = new Complex(0.,0.5*scale*timeStep*V[j]);
FDSolution: a.setL(j,mjh); //Matrix elements
FDSolution: a.setD(j,pjp1.add(z));
FDSolution: a.setR(j,mjh);
FDSolution: b.setL(j,pjh); //Right hand side
FDSolution: b.setD(j,mjp1.sub(z));
FDSolution: b.setR(j,pjh);
FDSolution: }
FDSolution:
FDSolution: c=b.dot(h); //Right hand side with
FDSolution: c[0]=c[0].add(b.getL(0).mul(h[n])); // with periodicity
FDSolution: c[n]=c[n].add(b.getR(n).mul(h[0]));
FDSolution:
FDSolution: hp=a.solve3(c); //Solve linear problem
FDSolution: for (int j=0; j<=n; j++){ //Plot norm & real part
FDSolution: fp[j]=hp[j].norm();
FDSolution: gp[j]=hp[j].re();
FDSolution: }
FDSolution:
FDSolution: } else if(pde.equals(jbone.BSCHOLES)
FDSolution: && scheme.equals(jbone.EUNAIVE)){
FDSolution: //=======================================================================
FDSolution: // FD - Black-Scholes - European Vanilla Put option - not normalized
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution: double disperCo= runData.getParamValue(runData.dispersionNm);
FDSolution:
FDSolution: double strike = runData.getParamValue(runData.icPositionNm);
FDSolution: double sigmaSq = diffusCo; //Volatility square
FDSolution: double rate = velocity; //Interest rate
FDSolution: double divid = disperCo; //Dividend
FDSolution:
FDSolution: fp[0]=strike*Math.exp(-rate*time); //Boundary condition
FDSolution: for (int j=1; j<n; j++) { //Explicit 2 levels
FDSolution: fp[j]=f[j+1]* 0.5*timeStep*(sigmaSq*j*j +(rate-divid)*j)
FDSolution: +f[j ]*(1.0-timeStep*(sigmaSq*j*j + rate ))
FDSolution: +f[j-1]* 0.5*timeStep*(sigmaSq*j*j -(rate-divid)*j);
FDSolution: }
FDSolution: fp[n]=fp[n-1]; //Boundary condition
FDSolution:
FDSolution: } else if(pde.equals(jbone.BSCHOLES)
FDSolution: && scheme.equals(jbone.EUPUT)){
FDSolution: //=======================================================================
FDSolution: // FD - Black-Scholes - European Vanilla Put option - normalized
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FDSolution: double disperCo= runData.getParamValue(runData.dispersionNm);
FDSolution:
FDSolution: double strike = runData.getParamValue(runData.icPositionNm);
FDSolution: double sigmaSq = diffusCo; //Volatility square
FDSolution: double rate = velocity; //Interest rate
FDSolution: double divid = disperCo; //Dividend
FDSolution: int j;
FDSolution: double x0, x1, f0, f1, xxi; //Change variables
FDSolution: double tau = 0.5*sigmaSq*time; // f(x,t) ->
FDSolution: double dtau= 0.5*sigmaSq*timeStep; // fm(xx,tau)
FDSolution: double xx0 = Math.log(x[1]/strike); //Lognormal mesh
FDSolution: double dxx =(Math.log(x[n]/strike)-Math.log(x[1]/strike))/(n-1);
FDSolution: double k1 = 2*rate/sigmaSq;
FDSolution: double k2 = 2*(rate-divid)/sigmaSq;
FDSolution: double k2m1= k2-1.;
FDSolution: double k2p1= k2+1.;
FDSolution:
FDSolution: //===== Interpolate IC from (x,t) to lognormal variables (xx,tau)
FDSolution: if (time<=timeStep) {
FDSolution: if(true) { //Initialize fm[] as put-option
FDSolution: for (int i=1; i<=n; i++) {
FDSolution: xxi=xx0+(i-1)*dxx;
FDSolution: fm[i]=Math.max(0., Math.exp(0.5*k2m1*xxi) -
FDSolution: Math.exp(0.5*k2p1*xxi) );
FDSolution: }
FDSolution: } else { //Interpolate fm[] from IC in f[]
FDSolution: j=1; ; x0=xx0;
FDSolution: f0=f[1]/strike*Math.exp(0.5*k2m1*xx0);
FDSolution: x1=x0; f1=f0;
FDSolution: for (int i=1; i<n; i++) { //Loop over lognormal mesh index
FDSolution: xxi=xx0+(i-1)*dxx; // find x0,x1 | x[j] < xxi < x[j+1]
FDSolution: while (xxi>=x1) { j++; x0=x1; f0=f1; x1=Math.log(x[j]/strike); }
FDSolution: f1=f[j]/strike*Math.exp(0.5*k2m1*x1);
FDSolution: fm[i]= f0 + (xxi-x0)*(f1-f0)/(x1-x0);
FDSolution: }
FDSolution: fm[n]= fm[n-1] + dxx*(fm[n-1]-fm[n-2]);
FDSolution: }
FDSolution: } else { //Retrieve fm[] from previous time step
FDSolution: }
FDSolution:
FDSolution: //===== Solve diffusion equation an explicit 2 levels scheme
FDSolution: double D = dtau/(dxx*dxx);
FDSolution: for (int i=2; i<n; i++)
FDSolution: f[i]= fm[i] + D*(fm[i+1]-2.*fm[i]+fm[i-1]);
FDSolution: f[1]= Math.exp(0.5*k2m1*xx0+0.25*k2m1*k2m1*tau); //Boundary conditions
FDSolution: f[n]= f[n-1];
FDSolution:
FDSolution: //===== Interpolate back from lognormal to financial mesh variables
FDSolution: fp[0]=strike*Math.exp(-rate*time); //Analytically
FDSolution: j=1; x0=x[0]; x1=x0; f0=fp[0];
FDSolution: xxi=xx0; f1=f[1]*strike*Math.exp(-0.5*k2m1*xxi-(0.25*k2m1*k2m1+k1)*tau);
FDSolution: for (int i=1; i<n; i++) { //Loop over financial mesh index
FDSolution: while (x[i]>=x1){
FDSolution: j++;x0=x1;f0=f1;xxi=xx0+(j-1)*dxx;x1=strike*Math.exp(xxi);
FDSolution: }
FDSolution: f1=f[j]*strike*Math.exp(-0.5*k2m1*xxi-(0.25*k2m1*k2m1+k1)*tau);
FDSolution: fp[i]= f0 +(x[i]-x0)/(x1-x0)*(f1-f0); //Lin interpol in x
FDSolution: // double gg=-0.5*k2m1;double dx=xxi-Math.log(x[i]/strike);//Expand exp()
FDSolution: // gp[i]= f0*(1.+gg*dx +0.5*gg*gg*dx*dx); // is worse;
FDSolution: }
FDSolution: xxi=Math.log(x[n]/strike);
FDSolution: fp[n]=f[n]*strike*Math.exp(-0.5*k2m1*xxi-(0.25*k2m1*k2m1+k1)*tau);
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_01)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.01
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.01.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_02)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.02
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.02.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_03)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.03
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.03.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_04)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.04
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.04.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_05)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.05
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.05.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_06)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.06
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.06.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_07)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.07
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.07.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_08)){
FDSolution: //=======================================================================
FDSolution: // FD - Exercise 2.08
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.08.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.MY_SCHEME)){
FDSolution: //=======================================================================
FDSolution: // FD - Project 1
FDSolution: //=======================================================================
FDSolution: //INCLUDE_prj2.01.java_INCLUDE//
FDSolution:
FDSolution: } else if(pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.MY_SCHEME_2)){
FDSolution: //=======================================================================
FDSolution: // FD - Project 2
FDSolution: //=======================================================================
FDSolution: //INCLUDE_prj2.02.java_INCLUDE//
FDSolution:
FDSolution: } // if
FDSolution:
FDSolution: if(isDefined && isForward)
FDSolution: shiftLevels();
FDSolution: return isDefined;
FDSolution: } // next
FDSolution:
FDSolution:
FDSolution: /** Take the solution backward one step to initialize schemes
FDSolution: with 3 time levels.
FDSolution: The equation is set by the internal variable pde and
FDSolution: the numerical scheme by the internal variable scheme
FDSolution: @param runData List of run parameters
FDSolution: @see RunData
FDSolution: @return True if an initialisation scheme is implemented for that equation
FDSolution: ****************************************************************************/
FDSolution: public boolean previous(RunData runData, PhysData physData){
FDSolution: int n=f.length-1;
FDSolution: boolean isDefined=false;
FDSolution:
FDSolution: if(pde.equals(jbone.ADVECTION) && scheme.equals(jbone.EXP3LVL) ||
FDSolution: pde.equals(jbone.KDV) && scheme.equals(jbone.EXP3LVL) ){
FDSolution:
FDSolution: //=======================================================================
FDSolution: // FD init - 3 levels - Common to all such schemes
FDSolution: //=======================================================================
FDSolution: for (int j=0; j<=n; j++) //Note how backward is implemented:
FDSolution: { fm[j]=0.;} // first reset the oldest level,
FDSolution: isForward = false; // then calculate a pseudo-forward
FDSolution: this.next(runData, physData); // solution without shifting
FDSolution: isForward = true; // levels.
FDSolution:
FDSolution: for (int j=0; j<=n; j++) { // Divide by (-2) to gives a scheme
FDSolution: fm[j]=f[j]-0.5*fp[j]; // equivalent to two levels backward
FDSolution: }
FDSolution:
FDSolution: } else if (pde.equals(jbone.ADVECTION) && scheme.equals(jbone.LPFROG)||
FDSolution: pde.equals(jbone.WAVE) && scheme.equals(jbone.FDTDABS)){
FDSolution:
FDSolution: //=======================================================================
FDSolution: // FDTD init - Advection/Wave - Leap-Frog
FDSolution: //=======================================================================
FDSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FDSolution: double velocity= runData.getParamValue(runData.velocityNm);
FDSolution: double frac = runData.getParamValue(runData.physDataValueNm);
FDSolution: double beta;
FDSolution: if (scheme.equals(jbone.LPFROG)) {
FDSolution: beta =timeStep*velocity/(dx[0]);
FDSolution: } else {
FDSolution: beta =timeStep/(dx[0]); // Initialize in vacuum
FDSolution: frac =1.0; // Propagate to the right
FDSolution: }
FDSolution:
FDSolution: for (int j=0; j<=n-1; j++) // Initialize pulse running direction
FDSolution: gm[j] = 0.5*frac*(f[j+1]+f[j]); // frac in [-1;1] -> [left; right]
FDSolution: gm[n]=0.5*frac*(f[0]+f[n]); // frac = 0 -> propagates in both
FDSolution:
FDSolution: for (int j=0; j<=n-1; j++) { //time=0.5*time_step
FDSolution: g[j] = gm[j]-0.5*beta*(f[j+1]-f[j]); }
FDSolution: g[n]=gm[n]-0.5*beta*(f[0]-f[n]);;
FDSolution:
FDSolution: } else if (pde.equals(jbone.EXERCISE)
FDSolution: && scheme.equals(jbone.EXC2_04)){
FDSolution: //=======================================================================
FDSolution: // FD init - Exersice 2.04
FDSolution: //=======================================================================
FDSolution: //INCLUDE_sol2.04.init_INCLUDE//
FDSolution:
FDSolution: }
FDSolution: return isDefined;
FDSolution: } // previous
FDSolution:
FDSolution: /** Tells whether the solution implements a option
FDSolution: @param option The the option to implement
FDSolution: @return True if the option is implemented
FDSolution: @see jbone#pdeNames
FDSolution: @see jbone#schemeNames
FDSolution: */
FDSolution: public boolean hasOption(String option){
FDSolution: try {
FDSolution: if(option.equals(jbone.ADVECTION))
FDSolution: return true;
FDSolution: //else if (option.equals(jbone.WAVE))
FDSolution: // return true;
FDSolution: else if (option.equals(jbone.BURGER))
FDSolution: return true;
FDSolution: else if (option.equals(jbone.KDV))
FDSolution: return true;
FDSolution: else if (option.equals(jbone.SCHROED))
FDSolution: return true;
FDSolution: else if (option.equals(jbone.BSCHOLES))
FDSolution: return true;
FDSolution: else if (option.equals(jbone.EXERCISE))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.EXP2LVL))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.EXP3LVL))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.IMP2LVL))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.EXPLAX))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.IMPLAX))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.LAXFRIED))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.LPFROG))
FDSolution: return true;
FDSolution: //else if (pde.equals(jbone.WAVE) && option.equals(jbone.FDTDABS))
FDSolution: // return true;
FDSolution: else if (pde.equals(jbone.BURGER) && option.equals(jbone.EXP2LVL))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.KDV) && option.equals(jbone.EXP3LVL))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.SCHROED) && option.equals(jbone.IMP2LVL))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.BSCHOLES) && option.equals(jbone.EUNAIVE))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.BSCHOLES) && option.equals(jbone.EUPUT))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_01))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_02))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_03))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_04))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_05))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_06))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_07))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC2_08))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXCB_01))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME))
FDSolution: return true;
FDSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME_2))
FDSolution: return true;
FDSolution: return false;
FDSolution: } catch (NullPointerException e){
FDSolution: return false;
FDSolution: } // try
FDSolution: } // hasOption
FDSolution:
FDSolution: } // class FDSolution
FEMSolution: /* $Id: FEMSolution.java,v 1.26 2001/01/15 14:22:21 andrej Exp $ */
FEMSolution:
FEMSolution: /******************************************************************************
FEMSolution: * FEMSolution -- contains the finite elements solver.
FEMSolution: * @see FluidSolution
FEMSolution: * @see Solution
FEMSolution: ******************************************************************************/
FEMSolution: public class FEMSolution extends FluidSolution{
FEMSolution:
FEMSolution: /** Creates a FEMSolution object.
FEMSolution: @param runData The run time parameters
FEMSolution: */
FEMSolution: public FEMSolution(RunData runData){
FEMSolution: super(runData);
FEMSolution: } // FEMSolution
FEMSolution:
FEMSolution: /** Advance the solution forward one step in time.
FEMSolution: The equation is set by the internal variable pde and
FEMSolution: the numerical scheme by the internal variable scheme
FEMSolution: @param runData List of run parameters
FEMSolution: @see RunData
FEMSolution: @return True if a scheme is implemented for that equation
FEMSolution: ****************************************************************************/
FEMSolution: public boolean next(RunData runData, PhysData physData) {
FEMSolution: int n=f.length-1;
FEMSolution: boolean isDefined=true;
FEMSolution:
FEMSolution: for (int j=0; j<=n; j++) { gp[j]=g[j]; } //Plot initial condition
FEMSolution:
FEMSolution: if(pde.equals(jbone.ADVECTION)
FEMSolution: && scheme.equals(jbone.TUNFEM)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Advection/diffusion, tunable integration t=1. FD
FEMSolution: // t=1/3 linear FEM
FEMSolution: // t=0. constant FEM
FEMSolution: //=======================================================================
FEMSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FEMSolution: double velocity= runData.getParamValue(runData.velocityNm);
FEMSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FEMSolution: double theta = runData.getParamValue(runData.timeImplicitNm);
FEMSolution: double tune = runData.getParamValue(runData.tuneIntegrNm);
FEMSolution:
FEMSolution: BandMatrix a = new BandMatrix(3, f.length);
FEMSolution: BandMatrix b = new BandMatrix(3, f.length);
FEMSolution: double[] c = new double[f.length];
FEMSolution:
FEMSolution: double alpha = timeStep*diffusCo/(dx[0]*dx[0]);//Are only constant for
FEMSolution: double beta = timeStep*velocity/(dx[0] );// homogeneous problems
FEMSolution: double h = dx[0];
FEMSolution: double htm = h*(1-tune)/4;
FEMSolution: double htp = h*(1+tune)/4;
FEMSolution:
FEMSolution: for (int j=0; j<=n; j++) {
FEMSolution: a.setL(j, htm +h*(-0.5*beta -alpha)* theta );
FEMSolution: a.setD(j,2*(htp +h*( alpha)* theta ));
FEMSolution: a.setR(j, htm +h*( 0.5*beta -alpha)* theta );
FEMSolution: b.setL(j, htm +h*(-0.5*beta -alpha)*(theta-1) );
FEMSolution: b.setD(j,2*(htp +h*( alpha)*(theta-1)));
FEMSolution: b.setR(j, htm +h*( 0.5*beta -alpha)*(theta-1) );
FEMSolution: }
FEMSolution:
FEMSolution: c=b.dot(f); //Right hand side
FEMSolution: c[0]=c[0]+b.getL(0)*f[n]; // with periodicity
FEMSolution: c[n]=c[n]+b.getR(n)*f[0];
FEMSolution:
FEMSolution: fp=a.solve3(c); //Solve linear problem
FEMSolution:
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.BSCHOLES) &&
FEMSolution: (scheme.equals(jbone.EUNAIVE)||
FEMSolution: scheme.equals(jbone.AMNAIVE) )){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Black-Scholes - American Vanilla Put option - not normalized
FEMSolution: //=======================================================================
FEMSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FEMSolution: double velocity= runData.getParamValue(runData.velocityNm);
FEMSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FEMSolution: double theta = runData.getParamValue(runData.timeImplicitNm);
FEMSolution: double tune = runData.getParamValue(runData.tuneIntegrNm);
FEMSolution:
FEMSolution: double sigmaSq = diffusCo; //Volatility square
FEMSolution: double rate = velocity; //Interest rate
FEMSolution: double strike = runData.getParamValue(runData.icPositionNm);
FEMSolution: double divid = runData.getParamValue(runData.dispersionNm);
FEMSolution:
FEMSolution: BandMatrix a = new BandMatrix(3, f.length); //Linear problem a*x=c
FEMSolution: BandMatrix b = new BandMatrix(3, f.length);
FEMSolution: double[] c = new double[f.length];
FEMSolution: double alp = timeStep*sigmaSq/2.; //PDE coefficients
FEMSolution: double bet =-timeStep*(rate-divid-sigmaSq);
FEMSolution: double gam = timeStep*rate;
FEMSolution: double h = dx[0];
FEMSolution: double t0m,t0,t0p,t1m,t1,t1p,t2m,t2,t2p; //Tunable quadra coeff
FEMSolution: t0p= h*0.25*(1-tune);
FEMSolution: t0m= h*0.25*(1-tune); t0=h*0.5*(1+tune);
FEMSolution: for (int j=0; j<=n; j++) {
FEMSolution: t1m=h*0.25*( -2*j-tune+1); t1=h*0.5*(tune-1);
FEMSolution: t1p=h*0.25*( 2*j-tune+1);
FEMSolution: t2m=h*0.25*(-4*j*j+4*j-tune-1); t2=h*0.5*(4*j*j+tune+1);
FEMSolution: t2p=h*0.25*(-4*j*j-4*j-tune-1);
FEMSolution: a.setL(j,t0m*(1.+gam* theta) +(t1m*bet +t2m*alp)* theta );
FEMSolution: a.setD(j,t0 *(1.+gam* theta) +(t1 *bet +t2 *alp)* theta );
FEMSolution: a.setR(j,t0p*(1.+gam* theta) +(t1p*bet +t2p*alp)* theta );
FEMSolution: b.setL(j,t0m*(1.+gam*(theta-1)) +(t1m*bet +t2m*alp)*(theta-1));
FEMSolution: b.setD(j,t0 *(1.+gam*(theta-1)) +(t1 *bet +t2 *alp)*(theta-1));
FEMSolution: b.setR(j,t0p*(1.+gam*(theta-1)) +(t1p*bet +t2p*alp)*(theta-1));
FEMSolution: }
FEMSolution: c=b.dot(f);
FEMSolution:
FEMSolution: if (scheme.equals(jbone.EUNAIVE)) { //===== European =======
FEMSolution: a.setL(0, 0.);a.setD(0,1.);a.setR(0,0.); //Dirichlet in-the-money
FEMSolution: c[0]=strike*Math.exp(-rate*time);
FEMSolution: a.setL(n,-1.);a.setD(n,1.);a.setR(n,0.);c[n]=0.; //Neuman out-the-money
FEMSolution: fp=a.ssor3(c,f);
FEMSolution:
FEMSolution: } else if (scheme.equals(jbone.AMNAIVE)) { //===== American =======
FEMSolution: double[] min = new double[f.length];
FEMSolution: double[] max = new double[f.length];
FEMSolution: for (int j=0; j<=n; j++) {
FEMSolution: min[j]=Math.max(strike-j*dx[0],0);
FEMSolution: max[j]=Double.POSITIVE_INFINITY;
FEMSolution: }
FEMSolution: a.setL(0, 0.);a.setD(0,1.);a.setR(0,0.); //Dirichlet in-the-money
FEMSolution: c[0]=strike*Math.exp(-rate*time);
FEMSolution: a.setL(n,-1.);a.setD(n,1.);a.setR(n,0.);c[n]=0.; //Neuman out-the-money
FEMSolution: double precision = strike*Math.pow(10.,-6); //Relative to strike
FEMSolution: int maxIter = 30;
FEMSolution: double w = 1.2; //Relaxation parameter
FEMSolution: fp=a.ssor3(c,f,min,max,precision,w,maxIter); //Projected-SSOR
FEMSolution: }
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.BSCHOLES) &&
FEMSolution: (scheme.equals(jbone.EUPUT)||
FEMSolution: scheme.equals(jbone.AMPUT) )){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Black-Scholes - American Vanilla Put option - normalized
FEMSolution: //=======================================================================
FEMSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FEMSolution: double velocity= runData.getParamValue(runData.velocityNm);
FEMSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FEMSolution: double theta = runData.getParamValue(runData.timeImplicitNm);
FEMSolution: double tune = runData.getParamValue(runData.tuneIntegrNm);
FEMSolution:
FEMSolution: double sigmaSq = diffusCo; //Volatility square
FEMSolution: double rate = velocity; //Interest rate
FEMSolution: double strike = runData.getParamValue(runData.icPositionNm);
FEMSolution: double divid = runData.getParamValue(runData.dispersionNm);
FEMSolution:
FEMSolution: int j;
FEMSolution: double x0, x1, xxi, f1,f0=0.; //Change variables
FEMSolution: double tau = 0.5*sigmaSq*time; // f(x,t) ->
FEMSolution: double dtau= 0.5*sigmaSq*timeStep; // fm(xx,tau)
FEMSolution: double xx0 = Math.log(x[1]/strike); //Lognormal mesh
FEMSolution: double dxx =(Math.log(x[n]/strike)-Math.log(x[1]/strike))/(n-1);
FEMSolution: double k1 = 2*rate/sigmaSq;
FEMSolution: double k2 = 2*(rate-divid)/sigmaSq;
FEMSolution: double k2m1= k2-1.;
FEMSolution: double k2p1= k2+1.;
FEMSolution:
FEMSolution: //===== Interpolate IC from (x,t) to lognormal variables (xx,tau)
FEMSolution: if (time<=timeStep) { //Initialize fm[] as put-option
FEMSolution: for (int i=0; i<=n; i++) {
FEMSolution: xxi=xx0+(i-1)*dxx;
FEMSolution: fm[i]=Math.max(0., Math.exp(0.5*k2m1*xxi) -
FEMSolution: Math.exp(0.5*k2p1*xxi) );
FEMSolution: }
FEMSolution: } else { //Retrieve fm[] from previous time step
FEMSolution: }
FEMSolution:
FEMSolution: //===== Solve diffusion equation with tunable finite elements
FEMSolution: BandMatrix a = new BandMatrix(3, f.length); //Linear problem a*x=c
FEMSolution: BandMatrix b = new BandMatrix(3, f.length);
FEMSolution: double[] c = new double[f.length];
FEMSolution: double htm = dxx*(1-tune)/4; //Tunable quadra coeff
FEMSolution: double htp = dxx*(1+tune)/4;
FEMSolution: double halpha = dtau/dxx;
FEMSolution: for (int i=0; i<=n; i++) {
FEMSolution: a.setL(i, htm -halpha* theta );
FEMSolution: a.setD(i,2*(htp +halpha* theta ));
FEMSolution: a.setR(i, htm -halpha* theta );
FEMSolution: b.setL(i, htm -halpha*(theta-1) );
FEMSolution: b.setD(i,2*(htp +halpha*(theta-1)));
FEMSolution: b.setR(i, htm -halpha*(theta-1) );
FEMSolution: }
FEMSolution: c=b.dot(fm);
FEMSolution: a.setL(0,0.);a.setD(0,1.);a.setR(0,0.);c[0]=0; //First equation idle
FEMSolution:
FEMSolution: if (scheme.equals(jbone.EUPUT)) { //===== European =======
FEMSolution: a.setL(1, 0.);a.setD(1,1.);a.setR(1,0.); //Dirichlet in-the-money
FEMSolution: c[1]=Math.exp(0.5*k2m1*xx0+0.25*k2m1*k2m1*tau);
FEMSolution: a.setL(n,-1.);a.setD(n,1.);a.setR(n,0.);c[n]=0.; //Neuman out-the-money
FEMSolution: f=a.ssor3(c,fm);
FEMSolution: f0=strike*Math.exp(-rate*time);
FEMSolution:
FEMSolution: } else if (scheme.equals(jbone.AMPUT)) { //===== American =======
FEMSolution: double[] min = new double[f.length];
FEMSolution: double[] max = new double[f.length];
FEMSolution: for (int i=0; i<=n; i++) {
FEMSolution: xxi=xx0+(i-1)*dxx;
FEMSolution: min[i]=Math.exp((0.25*k2m1*k2m1 +k1)*tau) *
FEMSolution: Math.max(0., Math.exp(0.5*k2m1*xxi) -
FEMSolution: Math.exp(0.5*k2p1*xxi) );
FEMSolution: max[i]=Double.POSITIVE_INFINITY;
FEMSolution: fm[i]=Math.max(min[i],f[i]); //IC interpolation error
FEMSolution: }
FEMSolution: a.setL(1, 0.);a.setD(1,1.);a.setR(1,0.); //Dirichlet in-the-money
FEMSolution: c[1]=Math.max(min[1],Math.exp(0.5*k2m1*xx0+0.25*k2m1*k2m1*tau));
FEMSolution: a.setL(n,-1.);a.setD(n,1.);a.setR(n,0.);c[n]=0.; //Neuman out-the-money
FEMSolution: double precision = strike*Math.pow(10.,-6); //Relative to strike
FEMSolution: int maxIter = 30;
FEMSolution: double w = 1.2; //Relaxation parameter
FEMSolution: f=a.ssor3(c,fm,min,max,precision,w,maxIter); //Projected-SSOR
FEMSolution: f0=strike;
FEMSolution: }
FEMSolution:
FEMSolution: //===== Interpolate back from lognormal to financial mesh variables
FEMSolution: j=1; x0=x[0]; x1=x0; fp[0]=f0;
FEMSolution: xxi=xx0; f1=f[1]*strike*Math.exp(-0.5*k2m1*xxi-(0.25*k2m1*k2m1+k1)*tau);
FEMSolution: for (int i=1; i<n; i++) { //Loop over financial mesh index
FEMSolution: while (x[i]>=x1){
FEMSolution: j++;x0=x1;f0=f1;xxi=xx0+(j-1)*dxx;x1=strike*Math.exp(xxi);
FEMSolution: }
FEMSolution: f1=f[j]*strike*Math.exp(-0.5*k2m1*xxi-(0.25*k2m1*k2m1+k1)*tau);
FEMSolution: fp[i]= f0 + (x[i]-x0)/(x1-x0)*(f1-f0);
FEMSolution: }
FEMSolution: xxi=Math.log(x[n]/strike);
FEMSolution: fp[n]=f[n]*strike*Math.exp(-0.5*k2m1*xxi-(0.25*k2m1*k2m1+k1)*tau);
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.EXC3_01)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Exercise 3.01
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_sol3.01.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.EXC3_02)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Exercise 3.02
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_sol3.02.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.EXC3_03)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Exercise 3.03
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_sol3.03.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.EXC3_04)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Exercise 3.04
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_sol3.04.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.EXC3_05)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Exercise 3.05
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_sol3.05.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.EXC3_07)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Exercise 3.07
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_sol3.07.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.EXCC_01)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Exercise C.01
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_solC.01.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.MY_SCHEME)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Project 1
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_prj3.01.java_INCLUDE//
FEMSolution:
FEMSolution: } else if(pde.equals(jbone.EXERCISE)
FEMSolution: && scheme.equals(jbone.MY_SCHEME_2)){
FEMSolution: //=======================================================================
FEMSolution: // FEM - Project 2
FEMSolution: //=======================================================================
FEMSolution: //INCLUDE_prj3.02.java_INCLUDE//
FEMSolution:
FEMSolution: } // if
FEMSolution:
FEMSolution: if(isDefined)
FEMSolution: shiftLevels();
FEMSolution: return isDefined;
FEMSolution: } // next
FEMSolution:
FEMSolution:
FEMSolution: /** Take the solution backward one step to initialize schemes
FEMSolution: with 3 time levels; not really appropriate in this context.
FEMSolution: @param runData List of run parameters
FEMSolution: @return False since it is not used here
FEMSolution: ****************************************************************************/
FEMSolution: public boolean previous(RunData runData, PhysData physData){ return false; }
FEMSolution:
FEMSolution: /** Tells whether the solution implements a option
FEMSolution: @param option The the option to implement
FEMSolution: @return True if the option is implemented
FEMSolution: @see jbone#pdeNames
FEMSolution: @see jbone#schemeNames
FEMSolution: */
FEMSolution: public boolean hasOption(String option){
FEMSolution: try {
FEMSolution: if(option.equals(jbone.ADVECTION))
FEMSolution: return true;
FEMSolution: else if(option.equals(jbone.BSCHOLES))
FEMSolution: return true;
FEMSolution: else if (option.equals(jbone.EXERCISE))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.ADVECTION) && option.equals(jbone.TUNFEM))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.BSCHOLES) && option.equals(jbone.AMPUT))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.BSCHOLES) && option.equals(jbone.EUPUT))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.BSCHOLES) && option.equals(jbone.AMNAIVE))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.BSCHOLES) && option.equals(jbone.EUNAIVE))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC3_01))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC3_02))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC3_03))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC3_04))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC3_05))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC3_07))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXCC_01))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME))
FEMSolution: return true;
FEMSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME_2))
FEMSolution: return true;
FEMSolution: return false;
FEMSolution: } catch (NullPointerException e){
FEMSolution: return false;
FEMSolution: } // try
FEMSolution: } // hasOption
FEMSolution:
FEMSolution: } // class FEMSolution
FFT: /* $Id: FFT.java,v 1.6 1999/08/19 14:33:30 andrej Exp $ */
FFT:
FFT: /***************************************************************************
FFT: * FFT - Object oriented Fast Fourier Transfrom package.
FFT: * Adapted with minimal changes to routines from "Numerical Recipes",
FFT: * W.H.Press et al., Cambridge University Press (1986)
FFT: ***************************************************************************/
FFT: public class FFT {
FFT: /* Array dimension */ public int size;
FFT: /* Label configuration space */ final static int inXSpace =-1;
FFT: /* Label Fourier space */ final static int inKSpace = 1;
FFT: /* Label complex number */ final static int bothParts = 0;
FFT: /* Label real part or 1st transform */ final static int firstPart = 1;
FFT: /* Label imag. part or 2nd transform */ final static int secondPart= 2;
FFT:
FFT: /* Current space location */ protected int location;
FFT: /* Dataset stored internally */ protected Complex[] data;
FFT: /* True if two real datasets */ protected boolean isTwoReals;
FFT: /* NR definition */ protected final int toFourier= +1;
FFT: /* NR definition */ protected final int toConfig = -1;
FFT:
FFT:
FFT: /** FFT object constructed from ONE complex array of dimension power(n,2).
FFT: @param arg Complex array
FFT: @param location Either inXSpace or inKSpace
FFT: ****************************************************************************/
FFT: public FFT(Complex arg[], int location) {
FFT: this.location = location; isTwoReals = false;
FFT: size = arg.length;
FFT: data = new Complex[size];
FFT: for (int k=0; k<size; k++) {
FFT: data[k] = arg[k];
FFT: }
FFT: }
FFT:
FFT:
FFT: /** FFT object constructed from ONE real array of dimension power(n,2).
FFT: @param arg Real array
FFT: @param location Either inXSpace or inKSpace
FFT: ****************************************************************************/
FFT: public FFT(double arg[], int location) {
FFT: this.location = location; isTwoReals = true;
FFT: size = arg.length;
FFT: data = new Complex[size];
FFT: for (int k=0; k<size; k++) {
FFT: data[k] = new Complex(arg[k], 0.);
FFT: }
FFT: }
FFT:
FFT:
FFT: /** FFT object constructed out of TWO complex arrays of dimension power(n,2).
FFT: @param arg1 First dataset
FFT: @param arg2 Second dataset
FFT: @param location Either inXSpace or inKSpace
FFT: ****************************************************************************/
FFT: public FFT(Complex arg1[], Complex arg2[], int location) {
FFT: this.location = location; isTwoReals = true;
FFT: size = arg1.length;
FFT: data = new Complex[size];
FFT: for (int k=0; k<size; k++) {
FFT: data[k] = arg1[k].add((Complex.i).mul(arg2[k]));
FFT: }
FFT: }
FFT:
FFT:
FFT: /** FFT object constructed out of TWO real arrays of dimension power(n,2).
FFT: @param arg1 First dataset
FFT: @param arg2 Second dataset
FFT: @param location Either inXSpace or inKSpace
FFT: ****************************************************************************/
FFT: public FFT(double arg1[], double arg2[], int location) {
FFT: this.location = location; isTwoReals = true;
FFT: size = arg1.length;
FFT: data = new Complex[size];
FFT: for (int k=0; k<size; k++) {
FFT: data[k] = new Complex(arg1[k],arg2[k]);
FFT: }
FFT: }
FFT:
FFT:
FFT: /** FFT object constructed from another FFT object
FFT: @param fft An FFT object to be duplicated
FFT: ****************************************************************************/
FFT: public FFT(FFT fft) {
FFT: this.size = fft.size;
FFT: this.location = fft.location;
FFT: this.isTwoReals = fft.isTwoReals;
FFT: for (int k=0; k<size; k++)
FFT: this.data[k]=fft.data[k];
FFT: }
FFT:
FFT:
FFT: /** Transform array to Fourier space.
FFT: @param index bothParts = complex number,
FFT: firstPart = real part or first transform if there is two,
FFT: secondPart = imaginary part or second transform.
FFT: @param units to scale with for physical results
FFT: @return Complex array in Fourier space.
FFT: If isTwoReals, the real and imaginary parts contain the result
FFT: of transforms of two real numbers.
FFT: @see FFT
FFT: ****************************************************************************/
FFT: public Complex[] getFromKSpace(int index, double units) {
FFT: Complex[] sol = new Complex[size];
FFT: if (location == inXSpace) {
FFT: if (isTwoReals) {
FFT: if (index == firstPart) {
FFT: Complex[] trash = new Complex[size];
FFT: twofft(data, sol, trash, size); }
FFT: else if (index == secondPart) {
FFT: Complex[] trash = new Complex[size];
FFT: twofft(data, trash, sol, size); }
FFT: else { //index == bothParts
FFT: for (int k=0; k<size; k++) sol[k]=data[k];
FFT: four1(sol, size, toFourier);
FFT: }
FFT: location = inXSpace; //Don't want to store two spectra inKSpace
FFT: for (int k=0; k<size; k++)
FFT: sol[k]=sol[k].scale(1./units);
FFT: } else { //
FFT: four1(data, size, toFourier);
FFT: location = inKSpace;
FFT: for (int k=0; k<size; k++)
FFT: sol[k]=data[k].scale(1./units);
FFT: }
FFT: } else {
FFT: for (int k=0; k<size; k++) sol[k]=data[k];
FFT: }
FFT: return sol;
FFT: }
FFT:
FFT:
FFT: /** Transform array to Fourier space and returns a real array depending
FFT: on the argument.
FFT: @param index firstPart = real part or first transform if there is two,
FFT: secondPart = imaginary part or second transform.
FFT: @param units to scale with for physical results
FFT: @return Real array in Fourier space.
FFT: @see FFT
FFT: ****************************************************************************/
FFT: public double[] getFromKSpacePart(int index, double units) {
FFT: Complex[] in = new Complex[size];
FFT: double[] out = new double[size];
FFT: in = this.getFromKSpace(index, units);
FFT: for (int k=0; k<size; k++) {
FFT: if (index==firstPart) out[k] = in[k].re();
FFT: else out[k] = in[k].im();
FFT: }
FFT: return out;
FFT: }
FFT:
FFT:
FFT: /** Transform array to configuration space.
FFT: @param units to scale with for physical results
FFT: @return Complex array in X space.
FFT: If isTwoReals, the real and imaginary parts contain the
FFT: result of two real transforms.
FFT: @see FFT
FFT: ****************************************************************************/
FFT: public Complex[] getFromXSpace(double units) {
FFT: if (location == inKSpace) {
FFT: four1(data, size, toConfig); //Inverse transfrom
FFT: for (int k=0; k<size; k++) // and normalization
FFT: data[k]=data[k].scale(units/(double)(size));
FFT: location = inXSpace;
FFT: }
FFT: return data;
FFT: }
FFT:
FFT:
FFT: /** Transform array to configuration space.
FFT: @param index firstPart = real part or first transform,
FFT: secondPart = iminary part or second transform.
FFT: @param units to scale with for physical results
FFT: @return Real or imaginary part of a complex array in Xspace.
FFT: If isTwoReals, the real and imaginary parts contain the result
FFT: of two real transforms.
FFT: @see FFT
FFT: ****************************************************************************/
FFT: public double[] getFromXSpacePart(int index, double units) {
FFT: Complex[] in = new Complex[size];
FFT: double[] out = new double[size];
FFT: in = this.getFromXSpace(units);
FFT: for (int k=0; k<size; k++) {
FFT: if (index==1) out[k] = in[k].re();
FFT: else out[k] = in[k].im();
FFT: }
FFT: return out;
FFT: }
FFT:
FFT:
FFT: /** Computes the convolution of two REAL numbers stored in the real and
FFT: imaginary part of the data[] array which has previously been initialized.
FFT: This is implement here WITHOUT any spectral expansion to show the
FFT: effect of ALIASING in nonlinear problems.
FFT: @param units to scale with for physical results
FFT: @return Result of the convolution (carefull ALIASED !)
FFT: @see FFT
FFT: ****************************************************************************/
FFT: public Complex[] aliasedConvolution(double units) {
FFT: Complex[] out = new Complex[size];
FFT: Complex[] re = new Complex[size];
FFT: Complex[] im = new Complex[size];
FFT:
FFT: four1(data, size, toFourier); //Transform 2 fields
FFT:
FFT: for (int k=0; k<size; k++) //Convolution =
FFT: data[k]=new Complex(data[k].re()*data[k].im()/units,0.);// Fourier product
FFT:
FFT: four1(data, size, toConfig); //Inverse transform
FFT:
FFT: for (int k=0; k<size; k++)
FFT: data[k]=data[k].scale(units/(double)(size)); //Rescale
FFT:
FFT: return data;
FFT: }
FFT:
FFT:
FFT:
FFT: /** Computes the convolution of two REAL numbers stored in the real and
FFT: imaginary part of the data[] array which has previously been initialized.
FFT: Here, the spectrum is temporarily extended to twice its original size
FFT: by padding the transforms with zeros.
FFT: @param units to scale with for physical results
FFT: @return Result of a correct convolution (with no aliasing).
FFT: @see FFT
FFT: ****************************************************************************/
FFT: public Complex[] expandedConvolution(double units) {
FFT: Complex[] dExp= new Complex[2*size];
FFT: Complex[] re= new Complex[2*size];
FFT: Complex[] im= new Complex[2*size];
FFT:
FFT: for (int k=0; k<2*size-1; k++) //Initialize
FFT: dExp[k]=new Complex(0.);
FFT:
FFT: dExp[0]=data[0];
FFT: for (int k=1; k<=size/2; k++) { //Expand array
FFT: dExp[k ]=data[k];
FFT: dExp[2*size-k]=data[size-k];
FFT: }
FFT:
FFT: four1(dExp, 2*size, toFourier); //Transform 2 fields
FFT:
FFT: for (int k=0; k<2*size; k++) //Convolution =
FFT: dExp[k]=new Complex(dExp[k].re()*dExp[k].im()/units,0.);// Fourier product
FFT:
FFT: four1(dExp, 2*size, toConfig); //Inverse transform
FFT:
FFT: data[0]=dExp[0];
FFT: for (int k=1; k<=size/2; k++) { //Rescale
FFT: data[k ]=dExp[k].scale(units/(double)(2*size));
FFT: data[size-k]=dExp[2*size-k].scale(1./(double)(2*size));
FFT: }
FFT: return data;
FFT: }
FFT:
FFT:
FFT: /** Fast Fourier transforms of two real arrays of dimension power(n,2).
FFT: Given a complex input array formed by packing two real arrays into
FFT: real and imaginary parts, this routine calls four1() and returns
FFT: two complex output arrays fft1[0:n-1] and fft2[0:n-1], each of complex
FFT: length n, which contain the discrete transforms of the respective
FFT: data arrays.
FFT: Adapted with minimal changes from "Numerical Recipes", W.H.Press et al.,
FFT: Cambridge University Press (1986)
FFT: ****************************************************************************/
FFT: protected void twofft(Complex datac[],
FFT: Complex fft1[], Complex fft2[], int n){
FFT: Complex h1 = new Complex();
FFT: Complex h2 = new Complex();
FFT: Complex c1 = new Complex(0.5, 0.0);
FFT: Complex c2 = new Complex(0.0,-0.5);
FFT:
FFT: for (int j=0;j<n;j++)
FFT: fft1[j] = datac[j]; //datac[] remains preserved
FFT:
FFT: four1(fft1,n,1); //Transform
FFT:
FFT: fft2[0]=new Complex(fft1[0].im(), 0.0);
FFT: fft1[0]=new Complex(fft1[0].re(), 0.0);
FFT: for (int j=1,k=n-j;j<=n/2;j++,k--) {
FFT: h1=fft1[j]; h1=h1.add(fft1[k].conj()); h1=h1.mul(c1);
FFT: h2=fft1[j]; h2=h2.sub(fft1[k].conj()); h2=h2.mul(c2);
FFT: fft1[j]=h1; fft1[k]=h1.conj();
FFT: fft2[j]=h2; fft2[k]=h2.conj();
FFT: }
FFT: }
FFT:
FFT:
FFT: /** Fast Fourier transform of one complex array of dimension power(n,2).
FFT: Replaces datac[0:n-1] by its discrete Fourier transform, if isign is 1,
FFT: or replaces datac[0:n-1] by n times its inverse discrete Fourier
FFT: transform, if isign is input as -1.
FFT: Adapted with minimal changes from "Numerical Recipes", W.H.Press et al.,
FFT: Cambridge University Press (1986)
FFT: ****************************************************************************/
FFT: protected void four1(Complex datac[], int nh, int isign){
FFT: int n,mmax,m,j,i,istep = 0;
FFT: double wtemp,wr,wpr,wpi,wi,theta;
FFT: double tempr,tempi;
FFT:
FFT: n = nh << 1; //Interface Complex to real
FFT: double[] data = new double[n+2];
FFT: for (int k=1;k<=nh;k++) {
FFT: data[2*k-1]=datac[k-1].re();
FFT: data[2*k ]=datac[k-1].im();
FFT: }
FFT:
FFT: j=1; //Bit reversal section
FFT: for (i=1;i<n;i+=2) {
FFT: if (j > i) { //Swap 2 complex numbers
FFT: tempr=data[j ];
FFT: tempi=data[j+1];
FFT: data[j ]=data[i ];
FFT: data[j+1]=data[i+1];
FFT: data[i ]=tempr;
FFT: data[i+1]=tempi;
FFT: }
FFT: m=n >> 1;
FFT: while (m >= 2 && j > m) {
FFT: j -= m;
FFT: m >>= 1;
FFT: }
FFT: j += m;
FFT: }
FFT: mmax=2; //Danielson-Lanczos section
FFT: while (n > mmax) { // executed log_2(n) times
FFT: istep=2*mmax;
FFT: theta=2*Math.PI/(isign*mmax); //Initialize trigonometric
FFT: wtemp=Math.sin(0.5*theta); // recurence
FFT: wpr = -2.0*wtemp*wtemp;
FFT: wpi=Math.sin(theta);
FFT: wr=1.0;
FFT: wi=0.0;
FFT: for (m=1;m<mmax;m+=2) { //Two nested inner loops
FFT: for (i=m;i<=n;i+=istep) {
FFT: j=i+mmax; //Danielson-Lanczos formula
FFT: tempr=wr*data[j ]-wi*data[j+1];
FFT: tempi=wr*data[j+1]+wi*data[j];
FFT: data[j ]=data[i ]-tempr;
FFT: data[j+1]=data[i+1]-tempi;
FFT: data[i ] += tempr;
FFT: data[i+1] += tempi;
FFT: }
FFT: wtemp=wr; //Trigonometric recurrence
FFT: wr=wtemp*wpr-wi*wpi+wr;
FFT: wi=wi*wpr+wtemp*wpi+wi;
FFT: }
FFT: mmax=istep;
FFT: }
FFT: for (int k=1;k<=nh;k++) //Back to Complex
FFT: datac[k-1] = new Complex(data[2*k-1], data[2*k]);
FFT: }
FFT:
FFT: }
FFTSolution: /* $Id: FFTSolution.java,v 1.30 2001/01/15 14:22:21 andrej Exp $ */
FFTSolution:
FFTSolution: /******************************************************************************
FFTSolution: * FFTsolution - Fast Fourier Transform solver
FFTSolution: * @see FluidSolution
FFTSolution: * @see Solution
FFTSolution: ******************************************************************************/
FFTSolution: public class FFTSolution extends FluidSolution{
FFTSolution: public FFTSolution(RunData runData){
FFTSolution: super(runData);
FFTSolution: } // FFTSolution
FFTSolution:
FFTSolution: /** FFT previous (or initial) step */ protected FFT keepFFT;
FFTSolution: /** FFT tool */ protected FFT toolFFT;
FFTSolution:
FFTSolution:
FFTSolution: /** Discretize the initial Shape function and initialize the moments
FFTSolution: @param function The initial shape to be approximated
FFTSolution: ****************************************************************************/
FFTSolution: public void discretize(ShapeFunction function){
FFTSolution: double boxLen = mesh.size()*mesh.interval(0);
FFTSolution: super.discretize(function);
FFTSolution: keepFFT = new FFT(f, FFT.inXSpace); //Load into FFT transformer
FFTSolution: } // discretize
FFTSolution:
FFTSolution:
FFTSolution: /** Advance the solution forward one step in time.
FFTSolution: The equation is set by the internal variable pde and
FFTSolution: the numerical scheme by the internal variable scheme
FFTSolution: @param runData List of run parameters
FFTSolution: @see RunData
FFTSolution: @return True if a scheme is implemented for that equation
FFTSolution: ****************************************************************************/
FFTSolution: public boolean next(RunData runData, PhysData physData) {
FFTSolution: boolean isForward = false;
FFTSolution: boolean isDefined = true;
FFTSolution:
FFTSolution: for (int i=0; i<f.length; i++) { gp[i]=g[i]; } // Plot initial value
FFTSolution:
FFTSolution: if(pde.equals(jbone.ADVECTION)
FFTSolution: && (scheme.equals(jbone.ALIASED)||
FFTSolution: scheme.equals(jbone.EXPAND ) )){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Advection/diffusion - No problem with aliasing if linear
FFTSolution: //=======================================================================
FFTSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FFTSolution: double velocity= runData.getParamValue(runData.velocityNm);
FFTSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FFTSolution: double disperCo= runData.getParamValue(runData.dispersionNm);
FFTSolution:
FFTSolution: int N = f.length; //Power of 2
FFTSolution: double boxLen = N*dx[0]; //Periodicity length
FFTSolution: Complex exp = new Complex();
FFTSolution: Complex total = new Complex();
FFTSolution: Complex[] s0 = new Complex[N];
FFTSolution:
FFTSolution: s0=keepFFT.getFromKSpace(FFT.firstPart,boxLen); //FFT real to KSpace
FFTSolution: // only once
FFTSolution: s[0] = s0[0];
FFTSolution: for (int m=1; m<N/2+1; m++) {
FFTSolution: double km = 2.*m*Math.PI/boxLen;
FFTSolution: Complex ikm1 = new Complex( 0. , km );
FFTSolution: Complex ikm2 = new Complex(-km*km, 0. );
FFTSolution: Complex ikm3 = new Complex( 0. ,-km*km*km);
FFTSolution: Complex advection = new Complex(ikm1.scale(velocity));
FFTSolution: Complex diffusion = new Complex(ikm2.scale(diffusCo));
FFTSolution: Complex dispersion= new Complex(ikm3.scale(disperCo));
FFTSolution: total=advection.add(diffusion);
FFTSolution: exp=(total.scale(time)).exp(); //Propagate directly
FFTSolution: s[m ] = s0[m].mul(exp); // from initial cond
FFTSolution: s[N-m] = s[m].conj(); //Here f is real
FFTSolution: }
FFTSolution: FFT ffts = new FFT(s,FFT.inKSpace); //Initialize Kspace
FFTSolution:
FFTSolution: f = ffts.getFromXSpacePart(FFT.firstPart,boxLen); //FFT back for plot
FFTSolution:
FFTSolution: } else if((pde.equals(jbone.KDV)||
FFTSolution: pde.equals(jbone.BURGER))
FFTSolution: && (scheme.equals(jbone.ALIASED)||
FFTSolution: scheme.equals(jbone.EXPAND ) )){
FFTSolution: //=======================================================================
FFTSolution: // FFT - KdV solitons and Burger's shocks - Aliasing here is an issue
FFTSolution: //=======================================================================
FFTSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
FFTSolution: double velocity= runData.getParamValue(runData.velocityNm);
FFTSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
FFTSolution: double disperCo= runData.getParamValue(runData.dispersionNm);
FFTSolution:
FFTSolution: int N = f.length; //Power of 2
FFTSolution: double boxLen = N*dx[0]; //Periodicity length
FFTSolution: Complex exp = new Complex();
FFTSolution: Complex nonlin = new Complex();
FFTSolution: Complex linear = new Complex();
FFTSolution: Complex total = new Complex();
FFTSolution: Complex[] s0 = new Complex[N];
FFTSolution:
FFTSolution: //----- Non-linear term: convolution
FFTSolution: s = keepFFT.getFromKSpace(FFT.bothParts,boxLen); //Current Spectrum
FFTSolution: toolFFT = new FFT(s,s,FFT.inKSpace); // for convolution
FFTSolution:
FFTSolution: if (scheme.equals(jbone.ALIASED)) //With-/out aliasing,
FFTSolution: sp = toolFFT.aliasedConvolution(boxLen); // use an FFT to
FFTSolution: else //scheme.equals(jbone.EXPAND) // calculate product,
FFTSolution: sp = toolFFT.expandedConvolution(boxLen); // FFT back to KSpace
FFTSolution:
FFTSolution:
FFTSolution: //----- Linear terms: complex terms in spectrum s
FFTSolution: s = keepFFT.getFromKSpace(FFT.bothParts,boxLen); //Current Spectrum
FFTSolution: linear= s[0];
FFTSolution: sp[0]=linear;
FFTSolution: for (int m=1; m<=N/2; m++) {
FFTSolution: double km = 2.*m*Math.PI/boxLen;
FFTSolution: Complex ikm1 = new Complex( 0. , km );
FFTSolution: Complex ikm2 = new Complex(-km*km, 0. );
FFTSolution: Complex ikm3 = new Complex( 0. ,-km*km*km);
FFTSolution: Complex advection = new Complex(ikm1.scale(velocity));
FFTSolution: Complex diffusion = new Complex(ikm2.scale(diffusCo));
FFTSolution: Complex dispersion= new Complex(ikm3.scale(disperCo));
FFTSolution: total=advection.add(diffusion);
FFTSolution: total=total.add(dispersion);
FFTSolution: exp=(total.scale(timeStep)).exp();
FFTSolution: linear = s[m].mul(exp);
FFTSolution: nonlin = sp[m].mul(ikm1.scale(0.5*timeStep));
FFTSolution: sp[m ] = linear.add(nonlin);
FFTSolution: if (m<N/2) sp[N-m] = sp[m].conj(); //For a real spectrum
FFTSolution: }
FFTSolution:
FFTSolution: keepFFT = new FFT(sp,FFT.inKSpace); //Spectrum is complete
FFTSolution: toolFFT = new FFT(sp,FFT.inKSpace); //inv FFT for plotting
FFTSolution: f=toolFFT.getFromXSpacePart(FFT.firstPart,boxLen);
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.EXC4_01)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Exercise 4.01
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_sol4.01.java_INCLUDE//
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.EXC4_02)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Exercise 4.02
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_sol4.02.java_INCLUDE//
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.EXC4_03)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Exercise 4.03
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_sol4.03.java_INCLUDE//
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.EXC4_04)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Exercise 4.04
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_sol4.04.java_INCLUDE//
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.EXC4_05)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Exercise 4.05
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_sol4.05.java_INCLUDE//
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.EXCD_01)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Exercise D.01
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_solD.01.java_INCLUDE//
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.MY_SCHEME)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Project 1
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_prj4.01.java_INCLUDE//
FFTSolution:
FFTSolution: } else if(pde.equals(jbone.EXERCISE)
FFTSolution: && scheme.equals(jbone.MY_SCHEME_2)){
FFTSolution: //=======================================================================
FFTSolution: // FFT - Project 2
FFTSolution: //=======================================================================
FFTSolution: //INCLUDE_prj4.02.java_INCLUDE//
FFTSolution:
FFTSolution: }
FFTSolution:
FFTSolution: if(isDefined && isForward)
FFTSolution: shiftLevels();
FFTSolution: return isDefined;
FFTSolution: } // next
FFTSolution:
FFTSolution:
FFTSolution: /** Starting procedure for FFT schemes.
FFTSolution: @param runData List of run parameters
FFTSolution: @see RunData
FFTSolution: @return True if an initialisation scheme is implemented for that equation
FFTSolution: ****************************************************************************/
FFTSolution: public boolean previous(RunData runData, PhysData physData){
FFTSolution: int n=f.length-1;
FFTSolution: boolean isDefined=false;
FFTSolution: return isDefined;
FFTSolution: } // previous
FFTSolution:
FFTSolution: /** Tells whether the solution implements a option
FFTSolution: @param option The the option to implement
FFTSolution: @return True if the option is implemented
FFTSolution: @see jbone#pdeNames
FFTSolution: @see jbone#schemeNames
FFTSolution: */
FFTSolution: public boolean hasOption(String option){
FFTSolution: try {
FFTSolution: if(option.equals(jbone.ADVECTION))
FFTSolution: return true;
FFTSolution: else if (option.equals(jbone.BURGER))
FFTSolution: return true;
FFTSolution: else if (option.equals(jbone.KDV))
FFTSolution: return true;
FFTSolution: else if (option.equals(jbone.EXERCISE))
FFTSolution: return true;
FFTSolution: else if ((pde.equals(jbone.ADVECTION)||
FFTSolution: pde.equals(jbone.KDV)||pde.equals(jbone.BURGER)) &&
FFTSolution: (option.equals(jbone.ALIASED)||option.equals(jbone.EXPAND)))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC4_01))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC4_02))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC4_03))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC4_04))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC4_05))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXCD_01))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME))
FFTSolution: return true;
FFTSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME_2))
FFTSolution: return true;
FFTSolution: return false;
FFTSolution: } catch (NullPointerException e){
FFTSolution: return false;
FFTSolution: } // try
FFTSolution: } // hasOption
FFTSolution:
FFTSolution: } // class FFTSolution
FluidSolution: /* $Id: FluidSolution.java,v 1.13 2000/06/16 07:03:57 johanh Exp $ */
FluidSolution:
FluidSolution: /*****************************************************************************
FluidSolution: * FluidSolution -- Is an abstract class at the top of all the solvers that
FluidSolution: * that operate directly on the fluid function.
FluidSolution: * @version $Revision: 1.23 $
FluidSolution: * @see Solution
FluidSolution: ******************************************************************************/
FluidSolution: abstract class FluidSolution extends Solution{
FluidSolution:
FluidSolution: /** True if needs to shiftLevels */
FluidSolution: protected boolean isForward = true;
FluidSolution:
FluidSolution: /** Creates a FluidSolution object.
FluidSolution: @param runData The run time parameters
FluidSolution: ****************************************************************************/
FluidSolution: public FluidSolution(RunData runData){
FluidSolution: super(runData);
FluidSolution: } // FluidSolution
FluidSolution:
FluidSolution: /** Discretize the initial Shape function and initialize the moments
FluidSolution: @param function The initial shape to be approximated
FluidSolution: ****************************************************************************/
FluidSolution: public void discretize(ShapeFunction function){
FluidSolution: for (int j = 0; j < mesh.size(); j++){ //Discretizing is very simple
FluidSolution: if (function.isComplex()){
FluidSolution: f[j] = function.getValue(mesh, j);
FluidSolution: h[j] = function.getValueC(mesh, j);
FluidSolution: f0[j]= h[j].re();
FluidSolution: } else {
FluidSolution: f[j] = function.getValue(mesh, j);
FluidSolution: f0[j] = f[j];
FluidSolution: }
FluidSolution: g[j] = 0.;
FluidSolution: }
FluidSolution: initialMoments[0]=calculateMoments(0); //Set conserved quantities
FluidSolution: initialMoments[1]=calculateMoments(1);
FluidSolution: rescale();
FluidSolution: } // discretize
FluidSolution:
FluidSolution: /** Internal function for copying new -> old
FluidSolution: ****************************************************************************/
FluidSolution: protected void shiftLevels() {
FluidSolution: for (int i=0;i<f.length;i++) {
FluidSolution: fm[i]=f[i]; dfm[i]=df[i]; f[i]=fp[i]; df[i]=dfp[i]; fp[i]=0.; dfp[i]=0.;
FluidSolution: gm[i]=g[i]; dgm[i]=dg[i]; g[i]=gp[i]; dg[i]=dgp[i]; gp[i]=0.; dgp[i]=0.;
FluidSolution: sm[i]=s[i]; hm[i]=h[i]; s[i]=sp[i]; h[i]=hp[i]; sp[i]=new Complex();
FluidSolution: hp[i]=new Complex();
FluidSolution: }
FluidSolution: }
FluidSolution:
FluidSolution: } // class FluidSolution
MCPDrawSolution: /* $Id: MCPDrawSolution.java,v 1.2 2000/06/16 07:03:57 johanh Exp $ */
MCPDrawSolution:
MCPDrawSolution: import java.awt.*;
MCPDrawSolution:
MCPDrawSolution: /** Solves the equation with a Monte Carlo method and draws the particles
MCPDrawSolution: @version $Revision: 1.2 $
MCPDrawSolution: @see Solution
MCPDrawSolution: @see MCPSolution
MCPDrawSolution: */
MCPDrawSolution: public class MCPDrawSolution extends MCPSolution {
MCPDrawSolution:
MCPDrawSolution: /** Creates a MCPSolution object. Note that discretize must be called
MCPDrawSolution: before next is called for the first time.
MCPDrawSolution: @param runData The run time parameters
MCPDrawSolution: @see Solution#next
MCPDrawSolution: @see Solution#discretize
MCPDrawSolution: */
MCPDrawSolution: public MCPDrawSolution(RunData runData){
MCPDrawSolution: super(runData);
MCPDrawSolution: } // MCPDrawSolution
MCPDrawSolution:
MCPDrawSolution: /** Plots the solution
MCPDrawSolution: @param plotArea The plot area
MCPDrawSolution: @param offScrImage The off screen image to draw on
MCPDrawSolution: @param headers Whether to draw headers
MCPDrawSolution: */
MCPDrawSolution: public void plot(Canvas plotArea, Image offScrImage, boolean headers){
MCPDrawSolution: if(!distributionUpToDate)
MCPDrawSolution: generateDistribution();
MCPDrawSolution: super.plot(plotArea, offScrImage, headers);
MCPDrawSolution:
MCPDrawSolution: int[] win=this.getWinSize();
MCPDrawSolution: xSize=win[0]; xOffset=win[2];
MCPDrawSolution: ySize=win[1]; xOffset=win[3];
MCPDrawSolution: Graphics g = offScrImage.getGraphics();
MCPDrawSolution: double dx = (xSize-2*xOffset)/(x_1-x_0);
MCPDrawSolution: double dy = (ySize-2*yOffset)/(y_1-y_0);
MCPDrawSolution: int[] x = new int[mesh.size() + 2];
MCPDrawSolution: int[] y = new int[mesh.size() + 2];
MCPDrawSolution: g.setColor(Color.red);
MCPDrawSolution: if (numberOfParticles>10) {
MCPDrawSolution: for(int i = 0; i < numberOfParticles; i++)
MCPDrawSolution: g.fillRect((int)Math.round((particlePosition[i]-x_0)*dx)+xOffset-1,
MCPDrawSolution: ySize -(i%(ySize-2*yOffset))-yOffset-2,3,3);
MCPDrawSolution: } else {
MCPDrawSolution: for(int i = 0; i < numberOfParticles; i++)
MCPDrawSolution: g.fillRect((int)Math.round((particlePosition[i]-x_0)*dx)+xOffset-1,
MCPDrawSolution: ySize -(i*8)-yOffset-12,8,8);
MCPDrawSolution: }
MCPDrawSolution: } // plot
MCPDrawSolution:
MCPDrawSolution: } // class MCPSolution
MCSDrawSolution: /* $Id: MCSDrawSolution.java,v 1.2 2000/06/16 07:03:57 johanh Exp $ */
MCSDrawSolution:
MCSDrawSolution: import java.awt.*;
MCSDrawSolution:
MCSDrawSolution: /** Solves the equation with a Monte Carlo sampling of possible configurations
MCSDrawSolution: and drawn on the plot area
MCSDrawSolution: @version $Revision: 1.34 $
MCSDrawSolution: @see Solution
MCSDrawSolution: @see MCSSolution
MCSDrawSolution: */
MCSDrawSolution: public class MCSDrawSolution extends MCSSolution {
MCSDrawSolution:
MCSDrawSolution: /** Creates a MCSSolution object. Note that discretize must be called
MCSDrawSolution: before next is called for the first time.
MCSDrawSolution: @param runData The run time parameters
MCSDrawSolution: @see Solution#next
MCSDrawSolution: @see Solution#discretize
MCSDrawSolution: */
MCSDrawSolution: public MCSDrawSolution(RunData runData){
MCSDrawSolution: super(runData);
MCSDrawSolution: } // MCSDrawSolution
MCSDrawSolution:
MCSDrawSolution: /** Plots the solution
MCSDrawSolution: @param plotArea The plot area
MCSDrawSolution: @param offScrImage The off screen image to draw on
MCSDrawSolution: @param headers Whether to draw headers
MCSDrawSolution: */
MCSDrawSolution: public void plot(Canvas plotArea, Image offScrImage, boolean headers){
MCSDrawSolution: if(!solutionUpToDate)
MCSDrawSolution: expectedRealisation();
MCSDrawSolution: super.plot(plotArea, offScrImage, headers);
MCSDrawSolution:
MCSDrawSolution: int[] win=this.getWinSize();
MCSDrawSolution: xSize=win[0]; xOffset=win[2];
MCSDrawSolution: ySize=win[1]; xOffset=win[3];
MCSDrawSolution: Graphics g = offScrImage.getGraphics();
MCSDrawSolution: double dx = (xSize-2*xOffset)/(x_1-x_0);
MCSDrawSolution: double dy = (ySize-2*yOffset)/(y_1-y_0);
MCSDrawSolution: int[] x = new int[mesh.size() + 2];
MCSDrawSolution: int[] y = new int[mesh.size() + 2];
MCSDrawSolution:
MCSDrawSolution: int dotX, dotY, dotSize;
MCSDrawSolution: g.setColor(Color.red);
MCSDrawSolution: for(int k=0; k<numberOfRealisations; k++) {
MCSDrawSolution: dotX=(int)Math.round((currentState[k]-x_0)*dx)+xOffset-1;
MCSDrawSolution: if (numberOfRealisations<10) {
MCSDrawSolution: dotSize=8; dotY=ySize-(k*dotSize)-yOffset-12;
MCSDrawSolution: } else {
MCSDrawSolution: dotSize=3; dotY=ySize-(k%(ySize-2*yOffset))-yOffset-2;
MCSDrawSolution: }
MCSDrawSolution: g.fillRect(dotX,dotY,dotSize,dotSize);
MCSDrawSolution: }
MCSDrawSolution: } // plot
MCSDrawSolution:
MCSDrawSolution: } // class MCSSolution
MCPSolution: /* $Id: MCPSolution.java,v 1.34 2001/01/15 14:22:22 andrej Exp $ */
MCPSolution:
MCPSolution: import java.awt.*;
MCPSolution: import java.util.Random;
MCPSolution:
MCPSolution:
MCPSolution:
MCPSolution: /** Solves the equation with a Monte Carlo method.
MCPSolution: @version $Revision: 1.34 $
MCPSolution: @see Solution
MCPSolution: */
MCPSolution: public class MCPSolution extends ParticleSolution{
MCPSolution:
MCPSolution: /** Creates a MCPSolution object. Note that discretize must be called
MCPSolution: before next is called for the first time.
MCPSolution: @param runData The run time parameters
MCPSolution: @see Solution#next
MCPSolution: @see Solution#discretize
MCPSolution: */
MCPSolution: public MCPSolution(RunData runData){
MCPSolution: super(runData);
MCPSolution: } // MCPSolution
MCPSolution:
MCPSolution:
MCPSolution: /** Advance the solution forward one step in time.
MCPSolution: The equation is set by the internal variable pde and
MCPSolution: the numerical scheme by the internal variable scheme
MCPSolution: @param runData List of run parameters
MCPSolution: @see RunData
MCPSolution: @return True if a scheme is implemented for that equation
MCPSolution: */
MCPSolution: public boolean next(RunData runData, PhysData physData) {
MCPSolution: distributionUpToDate = false;
MCPSolution: boolean isDefined = true;
MCPSolution: for (int i=0; i<f.length; i++) { gp[i]=g[i]; } //Plot initial condition
MCPSolution:
MCPSolution: double[] lim = {mesh.point(0),
MCPSolution: mesh.point(mesh.size()-1) + dx[0]};
MCPSolution:
MCPSolution: if(pde.equals(jbone.ADVECTION) &&
MCPSolution: (scheme.equals(jbone.MCNORM))){
MCPSolution: //=========================================================================
MCPSolution: // MCP - Advection/diffusion
MCPSolution: //=========================================================================
MCPSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
MCPSolution: double velocity= runData.getParamValue(runData.velocityNm);
MCPSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
MCPSolution:
MCPSolution: for(int j = 0; j < numberOfParticles; j++){
MCPSolution: particlePosition[j] += velocity * timeStep +
MCPSolution: random.nextGaussian() *
MCPSolution: Math.sqrt(2 * diffusCo * timeStep);
MCPSolution: // Periodic boundary condition; the last point is not a solution point
MCPSolution: while(particlePosition[j] > lim[1])
MCPSolution: particlePosition[j] -= lim[1] - lim[0];
MCPSolution: while(particlePosition[j] < lim[0])
MCPSolution: particlePosition[j] += lim[1] - lim[0];
MCPSolution: } // for
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXC5_01)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise 5.01
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_sol5.01.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXC5_02)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise 5.02
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_sol5.02.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXC5_03)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise 5.03
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_sol5.03.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXC5_04)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise 5.04
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_sol5.04.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXC5_05)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise 5.05
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_sol5.05.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXC5_06)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise 5.06
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_sol5.06.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXC5_07)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise 5.07
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_sol5.07.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(scheme.equals(jbone.EXCE_01)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Exercise E.01
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_solE.01.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(pde.equals(jbone.EXERCISE)
MCPSolution: && scheme.equals(jbone.MY_SCHEME)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Project 1
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_prj5.01.java_INCLUDE//
MCPSolution:
MCPSolution: } else if(pde.equals(jbone.EXERCISE)
MCPSolution: && scheme.equals(jbone.MY_SCHEME_2)){
MCPSolution: //=======================================================================
MCPSolution: // MCP - Project 2
MCPSolution: //=======================================================================
MCPSolution: //INCLUDE_prj5.02.java_INCLUDE//
MCPSolution:
MCPSolution: } // if
MCPSolution: return isDefined;
MCPSolution: } // next
MCPSolution:
MCPSolution: /** Tells whether the solution implements a option
MCPSolution: @param option The the option to implement
MCPSolution: @return True if the option is implemented
MCPSolution: @see jbone#pdeNames
MCPSolution: @see jbone#schemeNames
MCPSolution: */
MCPSolution: public boolean hasOption(String option){
MCPSolution: try {
MCPSolution: if (option.equals(jbone.ADVECTION))
MCPSolution: return true;
MCPSolution: else if (option.equals(jbone.EXERCISE))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.ADVECTION) &&
MCPSolution: option.equals(jbone.MCNORM) )
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_01))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_02))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_03))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_04))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_05))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_06))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_07))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXCE_01))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME))
MCPSolution: return true;
MCPSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME_2))
MCPSolution: return true;
MCPSolution: return false;
MCPSolution: } catch (NullPointerException e){
MCPSolution: return false;
MCPSolution: } // try
MCPSolution: } // hasOption
MCPSolution:
MCPSolution: /** Plots the solution
MCPSolution: @param plotArea The plot area
MCPSolution: @param offScrImage The off screen image to draw on
MCPSolution: @param headers Whether to draw headers
MCPSolution: */
MCPSolution: public void plot(Canvas plotArea, Image offScrImage, boolean headers){
MCPSolution: if(!distributionUpToDate)
MCPSolution: generateDistribution();
MCPSolution: super.plot(plotArea, offScrImage, headers);
MCPSolution: } // plot
MCPSolution:
MCPSolution: } // class MCPSolution
MCSSolution: /* $Id: MCSSolution.java,v 1.34 2001/01/15 14:22:22 andrej Exp $ */
MCSSolution:
MCSSolution: import java.awt.*;
MCSSolution: import java.util.Random;
MCSSolution:
MCSSolution:
MCSSolution:
MCSSolution: /** Solves the equation with a Monte Carlo sampling of possible configurations
MCSSolution: @version $Revision: 1.34 $
MCSSolution: @see Solution
MCSSolution: */
MCSSolution: public class MCSSolution extends SamplingSolution{
MCSSolution:
MCSSolution: /** Creates a MCSSolution object. Note that discretize must be called
MCSSolution: before next is called for the first time.
MCSSolution: @param runData The run time parameters
MCSSolution: @see Solution#next
MCSSolution: @see Solution#discretize
MCSSolution: */
MCSSolution: public MCSSolution(RunData runData){
MCSSolution: super(runData);
MCSSolution: } // MCSSolution
MCSSolution:
MCSSolution:
MCSSolution: /** Advance the solution forward one step in time.
MCSSolution: The equation is set by the internal variable pde and
MCSSolution: the numerical scheme by the internal variable scheme
MCSSolution: @param runData List of run parameters
MCSSolution: @see RunData
MCSSolution: @return True if a scheme is implemented for that equation
MCSSolution: */
MCSSolution: public boolean next(RunData runData, PhysData physData) {
MCSSolution: solutionUpToDate = false;
MCSSolution: boolean isDefined = true;
MCSSolution:
MCSSolution: for (int j=0; j<f.length; j++) { gp[j]=g[j]; } //Plot initial condition
MCSSolution:
MCSSolution: double[] lim = {mesh.point(0),
MCSSolution: mesh.point(mesh.size()-1) + dx[0]};
MCSSolution:
MCSSolution: if(pde.equals(jbone.BSCHOLES) &&
MCSSolution: scheme.equals(jbone.MCLOGN) ){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Black-Scholes - European Vanilla Put option - lognormal walk
MCSSolution: //=======================================================================
MCSSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
MCSSolution: double velocity= runData.getParamValue(runData.velocityNm);
MCSSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
MCSSolution: double strike = runData.getParamValue(runData.icPositionNm);
MCSSolution: double sigmaSq = diffusCo; //Volatility square
MCSSolution: double rate = velocity; //Interest rate
MCSSolution:
MCSSolution: for(int k=0; k<numberOfRealisations; k++)
MCSSolution: currentState[k]+= currentState[k]*( -rate*timeStep +
MCSSolution: random.nextGaussian()*sigmaSq*Math.sqrt(2*timeStep) );
MCSSolution:
MCSSolution: } else if(pde.equals(jbone.BSCHOLES) &&
MCSSolution: scheme.equals(jbone.MCNORM) ){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Black-Scholes - European Vanilla Put option - normal walk
MCSSolution: //=======================================================================
MCSSolution: double timeStep= runData.getParamValue(runData.timeStepNm);
MCSSolution: double velocity= runData.getParamValue(runData.velocityNm);
MCSSolution: double diffusCo= runData.getParamValue(runData.diffusionNm);
MCSSolution: double strike = runData.getParamValue(runData.icPositionNm);
MCSSolution: double sigmaSq = diffusCo; //Volatility square
MCSSolution: double rate = velocity; //Interest rate
MCSSolution:
MCSSolution: for(int k=0; k<numberOfRealisations; k++)
MCSSolution: currentState[k]+= strike*( -rate*timeStep +
MCSSolution: random.nextGaussian()*sigmaSq*Math.sqrt(2*timeStep) );
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXC5_01)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise 5.01
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_sol5.01.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXC5_02)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise 5.02
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_sol5.02.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXC5_03)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise 5.03
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_sol5.03.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXC5_04)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise 5.04
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_sol5.04.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXC5_05)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise 5.05
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_sol5.05.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXC5_06)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise 5.06
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_sol5.06.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXC5_07)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise 5.07
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_sol5.07.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(scheme.equals(jbone.EXCE_01)){
MCSSolution: //=======================================================================
MCSSolution: // MCS - Exercise E.01
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_solE.01.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(pde.equals(jbone.EXERCISE)
MCSSolution: && scheme.equals(jbone.MY_SCHEME)){
MCSSolution: //=======================================================================
MCSSolution: // MCP - Project 1
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_prj5.01.java_INCLUDE//
MCSSolution:
MCSSolution: } else if(pde.equals(jbone.EXERCISE)
MCSSolution: && scheme.equals(jbone.MY_SCHEME_2)){
MCSSolution: //=======================================================================
MCSSolution: // MCP - Project 2
MCSSolution: //=======================================================================
MCSSolution: //INCLUDE_prj5.02.java_INCLUDE//
MCSSolution:
MCSSolution: } // if
MCSSolution: return isDefined;
MCSSolution: } // next
MCSSolution:
MCSSolution: /** Tells whether the solution implements a option
MCSSolution: @param option The the option to implement
MCSSolution: @return True if the option is implemented
MCSSolution: @see jbone#pdeNames
MCSSolution: @see jbone#schemeNames
MCSSolution: */
MCSSolution: public boolean hasOption(String option){
MCSSolution: try {
MCSSolution: if (option.equals(jbone.BSCHOLES))
MCSSolution: return true;
MCSSolution: else if (option.equals(jbone.EXERCISE))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.BSCHOLES) &&
MCSSolution: (option.equals(jbone.MCNORM)|| option.equals(jbone.MCLOGN) ))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_01))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_02))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_03))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_04))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_05))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_06))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXC5_07))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.EXCE_01))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME))
MCSSolution: return true;
MCSSolution: else if (pde.equals(jbone.EXERCISE) && option.equals(jbone.MY_SCHEME_2))
MCSSolution: return true;
MCSSolution: return false;
MCSSolution: } catch (NullPointerException e){
MCSSolution: return false;
MCSSolution: } // try
MCSSolution: } // hasOption
MCSSolution:
MCSSolution: /** Plots the solution
MCSSolution: @param plotArea The plot area
MCSSolution: @param offScrImage The off screen image to draw on
MCSSolution: @param headers Whether to draw headers
MCSSolution: */
MCSSolution: public void plot(Canvas plotArea, Image offScrImage, boolean headers){
MCSSolution: if(!solutionUpToDate)
MCSSolution: expectedRealisation();
MCSSolution: super.plot(plotArea, offScrImage, headers);
MCSSolution: } // plot
MCSSolution:
MCSSolution: } // class MCSSolution
Mesh: /* $Id: Mesh.java,v 1.12 2000/06/06 12:58:57 andrej Exp $ */
Mesh:
Mesh: /******************************************************************************
Mesh: * Mesh - Discretization in one dimension
Mesh: ******************************************************************************/
Mesh: class Mesh {
Mesh:
Mesh: /** Coordinates */ double[] x;
Mesh: /** Interval sizes */ double[] dx;
Mesh:
Mesh: /** Constructor
Mesh: @param numberOfMeshPoints The number of mesh points
Mesh: @param meshLeft The left coordinate of the mesh
Mesh: @param meshLength The width of the mesh
Mesh: */
Mesh: public Mesh(int numberOfMeshPoints, double meshLeft, double meshLength) {
Mesh: double start = meshLeft;
Mesh: double end = meshLeft + meshLength;
Mesh:
Mesh: x = new double[numberOfMeshPoints]; //Targets
Mesh: dx = new double[numberOfMeshPoints];
Mesh:
Mesh: double dx0=(end-start)/((double)(numberOfMeshPoints)); //Rules
Mesh: for (int j=0; j<x.length; j++) {
Mesh: x[j]=start+dx0*(double)j;
Mesh: dx[j]=dx0;
Mesh: }
Mesh: }
Mesh:
Mesh: /** Grid Size
Mesh: @return Number of mesh points */
Mesh: public int size() { return x.length; }
Mesh:
Mesh: /** Coordinate values
Mesh: @return An Array with coordinates */
Mesh: public double[] points() { return x; }
Mesh:
Mesh: /** Coordinate value
Mesh: @param i The index of a coordinate
Mesh: @return The value of a coordinate */
Mesh: public double point(int i) { return x[i]; }
Mesh:
Mesh: /** Spacing between mesh points
Mesh: @return An Array with the interval sizes */
Mesh: public double[] intervals() { return dx; }
Mesh:
Mesh: /** Space to next grid point
Mesh: @param i The index of a coordinate
Mesh: @return The interval size */
Mesh: public double interval(int i) { return dx[i]; }
Mesh:
Mesh: /** Coordinate boundaries
Mesh: @return {min(x), max(x) */
Mesh: public double[] limits() {
Mesh: double[] lim = {x[0], x[0]};
Mesh: for (int i=1;i<x.length;i++) {
Mesh: if (x[i]<lim[0]) {lim[0]=x[i];};
Mesh: if (x[i]>lim[1]) {lim[1]=x[i];};
Mesh: }
Mesh: return lim;
Mesh: }
Mesh:
Mesh: /** For debugging */
Mesh: public String toString() {
Mesh: StringBuffer str = new StringBuffer();
Mesh: str.append("RegularMesh[").append(Double.toString(x[0]))
Mesh: .append(" , ").append(Double.toString(x[2]))
Mesh: .append(" .. ").append(Double.toString(x[x.length-1]))
Mesh: .append("]");
Mesh: return str.toString();
Mesh: }
Mesh:
Mesh: }//End of Mesh
MyChoice: /* $Id: MyChoice.java,v 1.10 2002/07/24 16:22:56 pde Exp $ */
MyChoice:
MyChoice: import java.awt.Choice;
MyChoice:
MyChoice: /** Class for handling the choices of the PDE type, scheme etc
MyChoice: @version $Revision: 1.10 $
MyChoice: */
MyChoice: class MyChoice extends Choice{
MyChoice:
MyChoice: /** List of the possible options of the choice */
MyChoice: String[] options;
MyChoice:
MyChoice: /** Constructor
MyChoice: @param item A list of the possible choices
MyChoice: */
MyChoice: public MyChoice(String[] item){
MyChoice: for(int i = 0; i < item.length; i++)
MyChoice: this.addItem(item[i]);
MyChoice: if (this.countItems()==0) this.addItem(" ");
MyChoice: select(0);
MyChoice: // Make a copy of the object
MyChoice: this.options = (String[])item.clone();
MyChoice: } // MyChoice
MyChoice:
MyChoice: /** Synchronize the list with the valid choices for a solution
MyChoice: @param solution The Solution to sync
MyChoice: */
MyChoice: public void sync(Solution solution){
MyChoice: String oldSelected = this.getSelectedItem();
MyChoice: // Clear the choices
MyChoice: // this.removeAll(); // Does not work; begin workaround
MyChoice: for(int i=0; i<this.options.length; i++) {
MyChoice: if(solution.hasOption(options[i])) {
MyChoice: this.insert(options[i],0);
MyChoice: i=10000;
MyChoice: }
MyChoice: }
MyChoice: for (int i=1; i<this.getItemCount(); i++) this.remove(i);
MyChoice: for (int i=1; i<this.getItemCount(); i++) this.remove(i);
MyChoice: for (int i=1; i<this.getItemCount(); i++) this.remove(i);
MyChoice: for (int i=1; i<this.getItemCount(); i++) this.remove(i);
MyChoice: long strt = System.currentTimeMillis(); //slow down selection
MyChoice: while (System.currentTimeMillis()-strt<200) { };
MyChoice: //End workaround
MyChoice:
MyChoice: boolean exist = false;
MyChoice: for(int i = 0; i < this.options.length; i++) {
MyChoice: if(solution.hasOption(options[i])) {
MyChoice: addItem(options[i]);
MyChoice: if (oldSelected.equals(options[i])) exist = true;
MyChoice: }
MyChoice: }
MyChoice: // if (this.countItems()>1) this.remove(0); //yields out-of-bound error
MyChoice: try { this.select(oldSelected); }
MyChoice: catch (IllegalArgumentException e) { this.select(1); }
MyChoice: } // sync
MyChoice: } // class MyChoice
MyChoice:
MyDoubleEditor: /* $Id: MyDoubleEditor.java,v 1.3 2000/04/11 16:08:40 johanh Exp $ */
MyDoubleEditor:
MyDoubleEditor: /** Editor of double values
MyDoubleEditor: @version $Revision: 1.3 $
MyDoubleEditor: */
MyDoubleEditor: public class MyDoubleEditor extends MyEditor{
MyDoubleEditor: /** Constructor
MyDoubleEditor: @param desc The description
MyDoubleEditor: @param value The current value
MyDoubleEditor: @param caller The calling class
MyDoubleEditor: */
MyDoubleEditor: public MyDoubleEditor(String desc, double value, MyEditorNotable caller){
MyDoubleEditor: super(desc, value, caller);
MyDoubleEditor: field.setText(String.valueOf(value));
MyDoubleEditor: } // MyDoubleEditor
MyDoubleEditor:
MyDoubleEditor: /** Parses the field and tryes to put in the value. If not
MyDoubleEditor: succeeded, this method will fail silently
MyDoubleEditor: */
MyDoubleEditor: protected void parseField() {
MyDoubleEditor: try {
MyDoubleEditor: value = Double.valueOf(field.getText()).doubleValue();
MyDoubleEditor: } catch (NumberFormatException e){
MyDoubleEditor: } // try
MyDoubleEditor: } // parseField
MyDoubleEditor:
MyDoubleEditor: } // class MyDoubleEditor
MyEditor: /* $Id: MyEditor.java,v 1.2 2000/04/11 16:07:58 johanh Exp $ */
MyEditor:
MyEditor: import java.awt.*;
MyEditor: import java.text.*;
MyEditor:
MyEditor: /** Class for editing run time parameters
MyEditor: @version $Revision: 1.2 $
MyEditor: */
MyEditor: abstract public class MyEditor extends Frame {
MyEditor: /** Description of the editor */
MyEditor: private String desc;
MyEditor:
MyEditor: /** The Set button */
MyEditor: private Button setButton;
MyEditor: /** The Cancel button */
MyEditor: private Button cancelButton;
MyEditor: /** The text field where the value is edited */
MyEditor: protected TextField field = new TextField(10);
MyEditor:
MyEditor: /** The value */
MyEditor: protected double value;
MyEditor:
MyEditor: /** The class calling */
MyEditor: private MyEditorNotable caller;
MyEditor:
MyEditor: /** Constructor
MyEditor: @param desc The description
MyEditor: @param value The current value
MyEditor: @param caller The calling class
MyEditor: */
MyEditor: public MyEditor(String desc, double value, MyEditorNotable caller){
MyEditor: this.desc = new String(desc);
MyEditor: this.value = value;
MyEditor: this.caller = caller;
MyEditor:
MyEditor: Panel p1 = new Panel();
MyEditor: p1.add(new Label(desc));
MyEditor: p1.add(field);
MyEditor: add("Center", p1);
MyEditor:
MyEditor: Panel p2 = new Panel();
MyEditor: p2.setLayout(new FlowLayout(FlowLayout.RIGHT));
MyEditor: cancelButton = new Button("Cancel");
MyEditor: setButton = new Button("Set");
MyEditor: p2.add(cancelButton);
MyEditor: p2.add(setButton);
MyEditor: add("South", p2);
MyEditor:
MyEditor: pack();
MyEditor: show();
MyEditor: } // MyEditor
MyEditor:
MyEditor: /** Read-only constructor for non-registered users
MyEditor: @param desc The description
MyEditor: @param value The current value
MyEditor: */
MyEditor: public MyEditor(String desc, double value){
MyEditor: this.desc = new String(desc);
MyEditor: this.value = value;
MyEditor:
MyEditor: Panel p1 = new Panel();
MyEditor: p1.add(new Label(desc));
MyEditor: p1.add(field);
MyEditor: add("Center", p1);
MyEditor:
MyEditor: Panel p2 = new Panel();
MyEditor: p2.setLayout(new FlowLayout(FlowLayout.RIGHT));
MyEditor: cancelButton = new Button("Cancel");
MyEditor: p2.add(new Label("Restricted: please login to edit"));
MyEditor: p2.add(cancelButton);
MyEditor: add("South", p2);
MyEditor:
MyEditor: pack();
MyEditor: show();
MyEditor: } // MyEditor
MyEditor:
MyEditor: /** Double-click on parameter requires editing
MyEditor: @deprecated
MyEditor: **************************************************************************/
MyEditor: public boolean action(Event event, Object arg) {
MyEditor: if (event.target == setButton || event.target == field) {
MyEditor: parseField();
MyEditor: caller.myEditorNotify(desc, value);
MyEditor: } // if
MyEditor: hide();
MyEditor: return true;
MyEditor: } // action
MyEditor:
MyEditor: /** Parses the field and tryes to put in the value. If not
MyEditor: succeeded, this method will fail silently
MyEditor: */
MyEditor: abstract protected void parseField();
MyEditor:
MyEditor: } // MyEditor
MyEditorNotable: /* $Id: MyEditorNotable.java,v 1.1 2000/04/10 13:10:06 johanh Exp $ */
MyEditorNotable:
MyEditorNotable: /** Interface suppying the methods needed for a calls to be
MyEditorNotable: notified by MyEditor
MyEditorNotable: @see MyEditor
MyEditorNotable: @version $Revision: 1.1 $
MyEditorNotable: */
MyEditorNotable: interface MyEditorNotable{
MyEditorNotable: /** Called by MyEditor
MyEditorNotable: @param desc The description of the property
MyEditorNotable: @param value The value to set
MyEditorNotable: @see MyEditor
MyEditorNotable: @see MyEditorNotable
MyEditorNotable: */
MyEditorNotable: public void myEditorNotify(String desc, double value);
MyEditorNotable: } // interface MyEditorNotable
MyNullEditor: /* $Id: MyNullEditor.java,v 1.10 2002/07/24 16:22:56 pde Exp $ */
MyNullEditor:
MyNullEditor: /** Editor demo preventing change of values
MyNullEditor: @version $Revision: 1.10 $
MyNullEditor: */
MyNullEditor: public class MyNullEditor extends MyEditor{
MyNullEditor: /** Constructor
MyNullEditor: @param desc The description
MyNullEditor: @param value The current value
MyNullEditor: */
MyNullEditor: public MyNullEditor(String desc, double value){
MyNullEditor: super(desc, value);
MyNullEditor: field.setText(String.valueOf(value));
MyNullEditor: } // MyNullEditor
MyNullEditor:
MyNullEditor: /** Idle */
MyNullEditor: protected void parseField() {
MyNullEditor: try {
MyNullEditor: // value = Double.valueOf(field.getText()).doubleValue();
MyNullEditor: } catch (NumberFormatException e){
MyNullEditor: } // try
MyNullEditor: }
MyNullEditor:
MyNullEditor: } // class MyNullEditor
MyIntEditor: /* $Id: MyIntEditor.java,v 1.4 2000/04/11 16:07:58 johanh Exp $ */
MyIntEditor:
MyIntEditor: /** Editor of int values
MyIntEditor: @version $Revision: 1.4 $
MyIntEditor: */
MyIntEditor: public class MyIntEditor extends MyEditor{
MyIntEditor: /** Constructor
MyIntEditor: @param desc The description
MyIntEditor: @param value The current value
MyIntEditor: @param caller The calling class
MyIntEditor: */
MyIntEditor: public MyIntEditor(String desc, double value, MyEditorNotable caller){
MyIntEditor: super(desc, value, caller);
MyIntEditor: field.setText(String.valueOf((int)value));
MyIntEditor: } // MyIntEditor
MyIntEditor:
MyIntEditor: /** Parses the field and tryes to put in the value. If not
MyIntEditor: succeeded, this method will fail silently
MyIntEditor: */
MyIntEditor: protected void parseField() {
MyIntEditor: try {
MyIntEditor: value = Integer.parseInt(field.getText());
MyIntEditor: } catch (NumberFormatException e){
MyIntEditor: } // try
MyIntEditor: } // parseField
MyIntEditor:
MyIntEditor: } // class MyIntEditor
ParticleSolution: /* $Id: ParticleSolution.java,v 1.23 2000/06/25 20:39:55 andrej Exp $ */
ParticleSolution:
ParticleSolution: import java.util.Random;
ParticleSolution:
ParticleSolution: /*****************************************************************************
ParticleSolution: * ParticleSolution -- Is an abstract class at the top of all the solvers that
ParticleSolution: * that push particles in phase space.
ParticleSolution: * @version $Revision: 1.23 $
ParticleSolution: * @see Solution
ParticleSolution: ******************************************************************************/
ParticleSolution: abstract class ParticleSolution extends Solution{
ParticleSolution:
ParticleSolution: /** The number of test particles in the simulation */
ParticleSolution: protected int numberOfParticles;
ParticleSolution: /** A vector with the position of the particles */
ParticleSolution: protected double[] particlePosition;
ParticleSolution: /** The amount of mass each test particle represents */
ParticleSolution: protected double particleDensity;
ParticleSolution: /** Whether the distribution function is up to date */
ParticleSolution: protected boolean distributionUpToDate;
ParticleSolution: /** Random number generator */
ParticleSolution: protected Random random;
ParticleSolution:
ParticleSolution: /** Creates a ParticleSolution object. The method discretize() must be called
ParticleSolution: before next() is called for the first time.
ParticleSolution: @param runData The run time parameters
ParticleSolution: @see Solution#next
ParticleSolution: @see Solution#discretize
ParticleSolution: ****************************************************************************/
ParticleSolution: public ParticleSolution(RunData runData){
ParticleSolution: super(runData);
ParticleSolution: numberOfParticles = runData.getParamValueInt(runData.walkersNm);
ParticleSolution: particlePosition = new double[numberOfParticles];
ParticleSolution: random = new Random();
ParticleSolution: } // ParticleSolution
ParticleSolution:
ParticleSolution: /** Take the solution backward one step to initialize 3 time level schemes;
ParticleSolution: not appropriate in this context.
ParticleSolution: @param runData List of run parameters
ParticleSolution: @return False since it is not used here
ParticleSolution: ****************************************************************************/
ParticleSolution: public boolean previous(RunData runData, PhysData physData){
ParticleSolution: return false;
ParticleSolution: } // previous
ParticleSolution:
ParticleSolution: /** Discretize the initial Shape function and initialize the moments
ParticleSolution: Initializes a new set of particles, the number being determined by
ParticleSolution: the numerical scheme parameter.
ParticleSolution: @param function The initial shape to be approximated
ParticleSolution: @see Solution#setScheme
ParticleSolution: ****************************************************************************/
ParticleSolution: public void discretize(ShapeFunction function){
ParticleSolution: int j;
ParticleSolution: for (j=0; j<mesh.size(); j++){ //Initial condition
ParticleSolution: f[j] = function.getValue(mesh, j);
ParticleSolution: f0[j] = f[j];
ParticleSolution: g[j] = 0.;
ParticleSolution: }
ParticleSolution: double[] fLim = super.limits(); //Min, max
ParticleSolution: double[] xLim = mesh.limits();
ParticleSolution: double meshDensity=(mesh.size()-1.)/(xLim[1]-xLim[0]);
ParticleSolution:
ParticleSolution: int jmax=0; //----- Detect Dirac function
ParticleSolution: for (j=0; j<mesh.size(); j++){
ParticleSolution: if (f[j]>f[jmax]) jmax=j;
ParticleSolution: }
ParticleSolution: if ( function.getValue(mesh,Math.max(jmax-1,0) ) < 0.1 *
ParticleSolution: function.getValue(mesh,jmax) &&
ParticleSolution: function.getValue(mesh,Math.min(jmax+1,mesh.size())) < 0.1 *
ParticleSolution: function.getValue(mesh,jmax) ){
ParticleSolution: initialMoments[0]=dx[jmax]*f[jmax];
ParticleSolution: initialMoments[1]=0.;
ParticleSolution: initialMoments[2]=0.;
ParticleSolution: particleDensity =
ParticleSolution: initialMoments[0]* meshDensity*meshDensity/numberOfParticles;
ParticleSolution:
ParticleSolution: for(j=0; j<mesh.size(); j++)
ParticleSolution: f[j]=(f[j]-fLim[0])/(fLim[1]-fLim[0]);
ParticleSolution:
ParticleSolution: for(j=0; j<numberOfParticles; j++)
ParticleSolution: particlePosition[j] = mesh.point(jmax);
ParticleSolution:
ParticleSolution: } else { //----- Regular function
ParticleSolution:
ParticleSolution: initialMoments[0]=calculateMoments(0); //Conserved quantities
ParticleSolution: initialMoments[1]=calculateMoments(1);
ParticleSolution: initialMoments[2]=calculateMoments(2);
ParticleSolution:
ParticleSolution: particleDensity =
ParticleSolution: initialMoments[0]* meshDensity*meshDensity/numberOfParticles;
ParticleSolution:
ParticleSolution: for(j=0; j<mesh.size(); j++) //Normalize to (0,1)
ParticleSolution: f[j]=(f[j]-fLim[0])/(fLim[1]-fLim[0]);
ParticleSolution:
ParticleSolution: double x;
ParticleSolution: distributionUpToDate = true;
ParticleSolution: for(j=0; j<numberOfParticles; j++){ // Randomize the positions
ParticleSolution: while(random.nextDouble() >
ParticleSolution: getValue(x = xLim[0]+(xLim[1]-xLim[0])*random.nextDouble())
ParticleSolution: );
ParticleSolution: particlePosition[j] = x;
ParticleSolution: } // for
ParticleSolution: }
ParticleSolution: distributionUpToDate = false;
ParticleSolution: rescale();
ParticleSolution: } // discretize
ParticleSolution:
ParticleSolution: /** Calculates the deviation from the m:th initial moment.
ParticleSolution: @param m The order of the moment
ParticleSolution: @return The deviation of the m:th moment from the beginning
ParticleSolution: ****************************************************************************/
ParticleSolution: public double momentsDeviation(int m){
ParticleSolution: if(!distributionUpToDate)
ParticleSolution: generateDistribution();
ParticleSolution: return super.momentsDeviation(m);
ParticleSolution: } // moments
ParticleSolution:
ParticleSolution: /** Calculates the limits of the solution.
ParticleSolution: @return A vector consisting of {min(solution), max(solution)}
ParticleSolution: ****************************************************************************/
ParticleSolution: public double[] limits(){
ParticleSolution: if(!distributionUpToDate)
ParticleSolution: generateDistribution();
ParticleSolution: return super.limits();
ParticleSolution: } // limits
ParticleSolution:
ParticleSolution: /** Gives the value of the function for an index
ParticleSolution: @param index The index for which to get the value
ParticleSolution: @return The value of the distribution function at a given index
ParticleSolution: ****************************************************************************/
ParticleSolution: public double getValue(int index) {
ParticleSolution: if(!distributionUpToDate)
ParticleSolution: generateDistribution();
ParticleSolution: return super.getValue(index);
ParticleSolution: } // getValue
ParticleSolution:
ParticleSolution: /** Linear interpolation of the solution. Assumes a uniform mesh.
ParticleSolution: @param arg Argument
ParticleSolution: @return A linear interpolation of the function for a given argument
ParticleSolution: ****************************************************************************/
ParticleSolution: public double getValue(double arg){
ParticleSolution: if(!distributionUpToDate)
ParticleSolution: generateDistribution();
ParticleSolution: return super.getValue(arg);
ParticleSolution: } // getValue
ParticleSolution:
ParticleSolution: /** Generates the distribution function from the set of test particles.
ParticleSolution: Assumes a uniform mesh.
ParticleSolution: ****************************************************************************/
ParticleSolution: protected void generateDistribution(){
ParticleSolution: int j;
ParticleSolution: for(j=0; j<mesh.size(); j++) f[j] = 0;
ParticleSolution: double[] xLim = mesh.limits();
ParticleSolution: for(j=0; j<numberOfParticles; j++){
ParticleSolution: // Lower mesh cell index
ParticleSolution: int index = (int)Math.floor((particlePosition[j] - xLim[0])
ParticleSolution: / (xLim[1] - xLim[0]) *
ParticleSolution: (mesh.size() - 1));
ParticleSolution:
ParticleSolution: if(index >= 0 && index <= mesh.size()){ //Particle is in plot region
ParticleSolution: if(index == mesh.size()) index--;
ParticleSolution: double x_low = mesh.point(index);
ParticleSolution: double x_high = mesh.point(index)+dx[0];
ParticleSolution: double a = (particlePosition[j]-x_low)/(x_high-x_low);
ParticleSolution: double cell_area_left = 1/dx[0]; // Assume uniform mesh
ParticleSolution: double cell_area_right = 1/dx[0];
ParticleSolution: f[index] += (1-a)* particleDensity/cell_area_left;
ParticleSolution: f[(index+1)%mesh.size()] += a * particleDensity/cell_area_right;
ParticleSolution: } // if
ParticleSolution: } // for
ParticleSolution: distributionUpToDate = true;
ParticleSolution: } // generateDistribution
ParticleSolution:
ParticleSolution: } // ParticleSolution
PhysData: /* $Id: PhysData.java,v 1.7 2001/01/15 14:22:22 andrej Exp $ */
PhysData:
PhysData: import java.awt.*;
PhysData: import java.text.*;
PhysData:
PhysData: /******************************************************************************
PhysData: * PhysData -- Is the class defining physical parameters (such as a potential)
PhysData: * @see Solution
PhysData: ******************************************************************************/
PhysData: public class PhysData{
PhysData:
PhysData: /** Potential */ public double[] potential;
PhysData:
PhysData: /** Creates a PhysData parameters object.
PhysData: @param runData The run time parameters
PhysData: */
PhysData: public PhysData(RunData runData){
PhysData: // CPU intensive initializations go here...
PhysData: } // Physics
PhysData:
PhysData: /** Evolves the physical data in time
PhysData: */
PhysData: public void update(double time){ // No evolution in the input data
PhysData: } // update
PhysData:
PhysData: /** Return the potential on the discrete mesh
PhysData: @return Array with discrete values of the potential
PhysData: */
PhysData: public double[] getPotential(RunData runData){
PhysData:
PhysData: Mesh mesh = runData.getMesh();
PhysData: double[] x = mesh.points();
PhysData: int n = mesh.size();
PhysData: int physDataCase = runData.getParamValueInt(runData.physDataCaseNm);
PhysData: double physDataValue = runData.getParamValue(runData.physDataValueNm);
PhysData:
PhysData: switch(physDataCase) {
PhysData:
PhysData: case 1: //Barrier = smooth Heavyside(x)
PhysData: //=======================================================================
PhysData: potential = new double[n]; //Create the potential array
PhysData: double tanhRange = 10.; //Set range of a in tanh(a)
PhysData: double tanhDx=2.*tanhRange/((double)n);
PhysData: double posBarrier=0.5*(x[n-1]-x[0]); //Barrier in the middle
PhysData: for(int k=0; k<n; k++){
PhysData: double a = (double)k*tanhDx-tanhRange;
PhysData: double tanh = (Math.exp(a)-Math.exp(-a))/(Math.exp(a)+Math.exp(-a));
PhysData: potential[k] = physDataValue*0.5*(1+tanh);
PhysData: }
PhysData: break;
PhysData:
PhysData: case 2: //Crystal = harmonic potential
PhysData: //=======================================================================
PhysData: potential = new double[n]; //Create the potential array
PhysData: for(int k=0; k<n; k++){
PhysData: double a = 2.*Math.PI*(double)k/(double)n;
PhysData: double cos = Math.cos(a);
PhysData: potential[k] = physDataValue*0.5*(1+cos);
PhysData: }
PhysData: break;
PhysData:
PhysData: case 3: //Your choice
PhysData: //=======================================================================
PhysData: //INCLUDE_sol2.6.phys_INCLUDE//
PhysData: break;
PhysData:
PhysData: default:
PhysData: break;
PhysData: }
PhysData: return potential;
PhysData: } // getPotential
PhysData:
PhysData: } // PhysData
PlotArea: /* $Id: PlotArea.java,v 1.4 2000/04/11 13:59:20 johanh Exp $ */
PlotArea:
PlotArea: import java.awt.*;
PlotArea:
PlotArea: /** Class used for the plot area
PlotArea: @version $Revision: 1.4 $
PlotArea: */
PlotArea: class PlotArea extends Canvas{
PlotArea:
PlotArea: /** The solution to plot */
PlotArea: private Solution solution;
PlotArea:
PlotArea: /** Whether to draw headers */
PlotArea: private boolean headers = true;
PlotArea:
PlotArea: /** Set the solution to plot
PlotArea: @param solution The solution to plot
PlotArea: */
PlotArea: public void setSolution(Solution solution){
PlotArea: this.solution = solution;
PlotArea: } // setSolution
PlotArea:
PlotArea: /** Toggles whether to draw headers
PlotArea: */
PlotArea: public void toggleHeaders() {
PlotArea: headers = !headers;
PlotArea: } // toggleHeaders
PlotArea:
PlotArea: /** Plot the graphics
PlotArea: @param g Graphics to plot
PlotArea: */
PlotArea: public void paint(Graphics g){
PlotArea: Dimension d;
PlotArea: try{
PlotArea: // Java 1.1
PlotArea: d = getSize();
PlotArea: } catch (NoSuchMethodError e) {
PlotArea: //Java1.0
PlotArea: d = size();
PlotArea: } // try
PlotArea: Image offScrImage = createImage(d.width, d.height);
PlotArea: solution.plot(this, offScrImage, headers);
PlotArea: g.drawImage(offScrImage, 0, 0, this);
PlotArea: } // paint
PlotArea:
PlotArea: /**
PlotArea: * No need to clear anything; just paint.
PlotArea: */
PlotArea: public void update(Graphics g) {
PlotArea: paint(g);
PlotArea: }
PlotArea:
PlotArea: } // PlotArea
RunData: /* $Id: RunData.java,v 1.26 2001/01/16 18:35:55 andrej Exp $ */
RunData:
RunData: import java.awt.*;
RunData: import java.text.*;
RunData: import java.applet.Applet;
RunData:
RunData: /** Class taking care of run time parameters
RunData: @version $Revision: 1.26 $
RunData: */
RunData: public class RunData extends List implements MyEditorNotable {
RunData:
RunData: /** Advection velocity */
RunData: public final static String velocityNm = "Velocity";
RunData: /** Diffusion constant */
RunData: public final static String diffusionNm = "Diffusion";
RunData: /** Dispersion constant*/
RunData: public final static String dispersionNm = "Dispersion";
RunData: /** Interval duration */
RunData: public final static String timeStepNm = "TimeStep";
RunData: /** Spatial resolution */
RunData: public final static String meshPointsNm = "MeshPoints";
RunData: /** Number of random walkers (particles, states) */
RunData: public final static String walkersNm = "Walkers";
RunData: /** Implicit time parameter*/
RunData: public final static String timeImplicitNm = "TimeImplicit";
RunData: /** FEM tunable integration*/
RunData: public final static String tuneIntegrNm = "TuneIntegr";
RunData: /** IC amplitude pulse1*/
RunData: public final static String icAmplitudeNm = "ICAmplitude";
RunData: /** IC position pulse1 */
RunData: public final static String icPositionNm = "ICPosition";
RunData: /** IC width pulse1 */
RunData: public final static String icWidthNm = "ICWidth";
RunData: /** IC wavelength pulse1 */
RunData: public final static String icWavelengthNm = "ICWavelength";
RunData: /** Final time to reach*/
RunData: public final static String runTimeNm = "RunTime";
RunData: /** Left coordinate of x-axis */
RunData: public final static String meshLeftNm = "MeshLeft";
RunData: /** Length of the domain */
RunData: public final static String meshLengthNm = "MeshLength";
RunData: /** Physical data switch */
RunData: public final static String physDataCaseNm = "PhysDataCase";
RunData: /** Physical data value */
RunData: public final static String physDataValueNm = "PhysDataValue";
RunData: /** Maximum length of the parameters strings */
RunData: private int maxLength;
RunData: /** The names of all the run time parameters */
RunData: public final static String[] paramNames = {
RunData: velocityNm, diffusionNm, dispersionNm, timeStepNm,
RunData: meshPointsNm, walkersNm, timeImplicitNm, tuneIntegrNm,
RunData: icAmplitudeNm, icPositionNm, icWidthNm, icWavelengthNm,
RunData: runTimeNm, meshLeftNm, meshLengthNm, physDataCaseNm,
RunData: physDataValueNm};
RunData: /** The type (integer or double) of the parameter */
RunData: private final static boolean[] paramIsInt = {
RunData: false, false, false, false,
RunData: true, true, false, false,
RunData: false, false, false, false,
RunData: false, false, false, true,
RunData: false};
RunData: /** The value of the parameter converted to double */
RunData: private double[] paramVal;
RunData: /** Display the parameters in the applet or not */
RunData: private boolean[] paramShow;
RunData:
RunData: /** Discretization mesh*/ Mesh mesh;
RunData: /** The calling class */ private RunDataNotable caller;
RunData: /** Dessert */ private boolean cookme = false;
RunData:
RunData: /** Creates the mesh from the parameters */
RunData: public void createMesh(){
RunData: int meshPoints = getParamValueInt(meshPointsNm);
RunData: double meshLeft = getParamValue(meshLeftNm);
RunData: double meshLength = getParamValue(meshLengthNm);
RunData: mesh = new Mesh(meshPoints, meshLeft, meshLength);
RunData: } // createMesh
RunData:
RunData: /** Constructor */
RunData: public RunData(RunDataNotable caller) {
RunData: super(paramNames.length, false); // for window height if not applet
RunData: this.caller = caller;
RunData: // initialize parameters
RunData: maxLength = 0; int i;
RunData: paramShow = new boolean[paramNames.length];
RunData: paramVal = new double[paramNames.length];
RunData: for(i=0; i<paramNames.length; i++) {
RunData: paramShow[i] = false;
RunData: if(paramNames[i].length() > maxLength)
RunData: maxLength = paramNames[i].length();
RunData: }
RunData: maxLength =+ 11;
RunData:
RunData: setParamValue(velocityNm,1.); // Advection velocity
RunData: setParamValue(diffusionNm,0.); // Diffusion constant
RunData: setParamValue(dispersionNm,0.); // Dispersion constan
RunData: setParamValue(timeStepNm,0.5); // Interval duration
RunData: setParamValue(meshPointsNm,64.); // Spatial resolution
RunData: setParamValue(walkersNm,100.); // Number of walkers
RunData: setParamValue(timeImplicitNm,0.7); // Implicit time para
RunData: setParamValue(tuneIntegrNm,0.333); // FEM tune integratio
RunData: setParamValue(runTimeNm,64.); // Final time to reac
RunData: setParamValue(icAmplitudeNm,1.); // IC amplitude pulse
RunData: setParamValue(icPositionNm,18.); // IC position pulse1
RunData: setParamValue(icWidthNm,8.); // IC width pulse1
RunData: setParamValue(icWavelengthNm,4.); // IC wavelength pulse1
RunData: setParamValue(meshLeftNm,0.); // Left coordinate
RunData: setParamValue(meshLengthNm,64.); // Length of domain
RunData: setParamValue(physDataCaseNm,1.); // Physical data switch
RunData: setParamValue(physDataValueNm,1.); // Physical data switch
RunData:
RunData: // update panel
RunData: this.setBackground(new Color(Integer.parseInt("88FFFF",16)));
RunData: setFont(new Font("Courier", Font.PLAIN, 12));
RunData: syncList();
RunData: createMesh();
RunData: } // RunData
RunData:
RunData: /** Returns index of a parameter
RunData: @param name String The parameter name
RunData: @return The parameter index */
RunData: final public int getParamIndex(String name){
RunData: for(int i=0; i<paramNames.length; i++)
RunData: if (name.equals(paramNames[i])) return i;
RunData: return -1;
RunData: }
RunData:
RunData: /** Returns the value a double parameter
RunData: @param name String The parameter name
RunData: @return The parameter value */
RunData: final public double getParamValue(String name){
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: if (name.equals(paramNames[i])) return paramVal[i];
RunData: }
RunData: // This should not happen
RunData: System.out.println("Unknown param '"+name+"' in RunData.getParamValue()!");
RunData: return 0.;
RunData: }
RunData:
RunData: /** Returns the value of an integer parameter
RunData: @param name String The parameter name
RunData: @return The parameter value */
RunData: final public int getParamValueInt(String name){
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: if (name.equals(paramNames[i])) {
RunData: return (int)paramVal[i];
RunData: }
RunData: }
RunData: // This should not happen
RunData: System.out.println("Unknown param '"+name+"' in RunData.getParamValue()!");
RunData: return 0;
RunData: }
RunData:
RunData: /** Sets the value of a parameter
RunData: @param name String The parameter name
RunData: @return The parameter value */
RunData: final public void setParamValue(String name, double value){
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: if (name.equals(paramNames[i])) {
RunData: paramVal[i]=value;
RunData: syncList();
RunData: return;
RunData: }
RunData: }
RunData: // This should not happen
RunData: System.out.println("Unknown param '"+name+"' in RunData.setParamValue()!");
RunData: return;
RunData: }
RunData:
RunData: /** Synchronizes the list with the current parameters */
RunData: final public void syncList(){
RunData: removeAll();
RunData: StringBuffer buffer;
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: if (paramShow[i]) {
RunData: buffer = new StringBuffer(maxLength);
RunData: buffer.append(paramNames[i]); pad(buffer, maxLength);
RunData: if (paramIsInt[i])
RunData: buffer.append(" = " + String.valueOf((int)paramVal[i]));
RunData: else
RunData: buffer.append(" = " + String.valueOf(paramVal[i]));
RunData: addItem(buffer.toString());
RunData: }
RunData: }
RunData: } // syncList
RunData:
RunData: /** Modify defaults parameters the HTML tags from the web page
RunData: @param a The applet to get the parameters from
RunData: */
RunData: final public void tagModify(Applet a) {
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: try { String tmp = a.getParameter(paramNames[i]);
RunData: if(tmp!=null) {
RunData: paramVal[i]=Double.valueOf(tmp).doubleValue();
RunData: paramShow[i]=true;
RunData: }
RunData: } catch (NumberFormatException e) {}
RunData: }
RunData: String tmp = a.getParameter("Cookme");
RunData: if(tmp!=null) cookme = false;
RunData: syncList();
RunData: } // TagModify
RunData:
RunData: /** Selects all the parameters for display in the applet
RunData: */
RunData: final public void displayAll() {
RunData: for(int i=0; i<paramNames.length; i++)
RunData: paramShow[i]=true;
RunData: syncList();
RunData: } // displayAll
RunData:
RunData: /** Selects only tag parameters for display in the applet
RunData: */
RunData: final public void displayTag(Applet a) {
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: if(a.getParameter(paramNames[i]) !=null)
RunData: paramShow[i]=true;
RunData: else
RunData: paramShow[i]=false;
RunData: }
RunData: syncList();
RunData: } // displayTag
RunData:
RunData: /** Copy the string padding it up with sapces to a specified length
RunData: @param d The output string
RunData: @param length The length to pad to
RunData: */
RunData: private void pad(StringBuffer d, int length){
RunData: while(d.length() < length)
RunData: d.append(" ");
RunData: } // pad
RunData:
RunData: /** Responds to the user actions through the mouse and buttons (Java1.0).
RunData: Yes, we know this sucks compare to Java 1.1, but we want to be
RunData: compatible with as many browsers as possible. There is a lot of
RunData: old stuff out there...
RunData: @deprecated
RunData: */
RunData: final public boolean action(Event e, Object arg) {
RunData: if(e.target instanceof RunData) {
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: if(stringPartCompare(paramNames[i], (String)arg)) {
RunData: if (!cookme) {
RunData: if (paramIsInt[i]) {
RunData: MyEditor ed=new MyIntEditor(paramNames[i],(int)paramVal[i],this);
RunData: } else {
RunData: MyEditor ed=new MyDoubleEditor(paramNames[i],paramVal[i],this);
RunData: }
RunData: } else { // cookme
RunData: if (paramIsInt[i]) {
RunData: MyEditor ed=new MyNullEditor(paramNames[i],(int)paramVal[i]);
RunData: } else {
RunData: MyEditor ed=new MyNullEditor(paramNames[i],paramVal[i]);
RunData: }
RunData: } // cookme
RunData: return true;
RunData: } // stringCompare
RunData: } // for
RunData: }
RunData: return false;
RunData: } // action
RunData:
RunData: /** Called by MyEditor
RunData: @param name The description of the property
RunData: @param value The value to set
RunData: @see MyEditor
RunData: @see MyEditorNotable
RunData: */
RunData: final public void myEditorNotify(String name, double value){
RunData: for(int i=0; i<paramNames.length; i++) {
RunData: if(paramNames[i].equals(name)) paramVal[i]=value;
RunData: }
RunData: syncList();
RunData: } // myEditorNote
RunData:
RunData: /** Internal method usedin action for camparing if two strings are
RunData: alike upon the length of the shortest string
RunData: @param s1 String 1 (the sortest)
RunData: @param s2 String 2
RunData: @return True if s1 == s2 upon the length of s1
RunData: */
RunData: final private boolean stringPartCompare(String s1, String s2) {
RunData: return s2.substring(0, s1.length()).equals(s1);
RunData: } // stringPartCompare
RunData:
RunData: /** Returns the mesh @return The mesh */
RunData: final public Mesh getMesh(){ return mesh; }
RunData:
RunData: } // class RunData
RunDataNotable: /* $Id: RunDataNotable.java,v 1.3 2000/06/25 20:39:55 andrej Exp $ */
RunDataNotable:
RunDataNotable: /** Interface suppying the methods needed for a calls to be
RunDataNotable: notified by RunData
RunDataNotable: @see MyEditor
RunDataNotable: @version $Revision: 1.3 $
RunDataNotable: */
RunDataNotable: interface RunDataNotable{
RunDataNotable: /** Invoked by RunData
RunDataNotable: @see RunData
RunDataNotable: */
RunDataNotable: public void runDataNotifyMesh();
RunDataNotable: /** Invoked by RunData
RunDataNotable: @see RunData
RunDataNotable: */
RunDataNotable: public void runDataNotifyWalkers();
RunDataNotable: } // interface RunDataWalkers
SamplingSolution: /* $Id: SamplingSolution.java,v 1.23 2000/06/25 20:39:55 andrej Exp $ */
SamplingSolution:
SamplingSolution: import java.util.Random;
SamplingSolution:
SamplingSolution: /*****************************************************************************
SamplingSolution: * SamplingSolution -- Is an abstract class at the top of all the solvers that
SamplingSolution: * that evolve an expected solution by sampling possible realisation.
SamplingSolution: * @version $Revision: 1.23 $
SamplingSolution: * @see Solution
SamplingSolution: ******************************************************************************/
SamplingSolution: abstract class SamplingSolution extends Solution{
SamplingSolution:
SamplingSolution: /** The number of independent realisations to evolve */
SamplingSolution: protected int numberOfRealisations;
SamplingSolution: /** The weight or probability associated with each initial state */
SamplingSolution: protected double realisationWeight;
SamplingSolution: /** A vector with the current state of each realisation */
SamplingSolution: protected double[] currentState;
SamplingSolution: /** Whether the expected solution is up to date */
SamplingSolution: protected boolean solutionUpToDate;
SamplingSolution: /** Random number generator */
SamplingSolution: protected Random random;
SamplingSolution:
SamplingSolution: /** Creates a SamplingSolution object. Note that discretize must be called
SamplingSolution: before next is called for the first time.
SamplingSolution: @param runData The run time parameters
SamplingSolution: @see Solution#next
SamplingSolution: @see Solution#discretize
SamplingSolution: ****************************************************************************/
SamplingSolution: public SamplingSolution(RunData runData){
SamplingSolution: super(runData);
SamplingSolution: numberOfRealisations = runData.getParamValueInt(runData.walkersNm);
SamplingSolution: currentState = new double[numberOfRealisations];
SamplingSolution: random = new Random();
SamplingSolution: } // SamplingSolution
SamplingSolution:
SamplingSolution: /** Take the solution backward one step to initialize schemes
SamplingSolution: with 3 time levels; not really appropriate in this context.
SamplingSolution: @param runData List of run parameters
SamplingSolution: @return False since it is not used here
SamplingSolution: ****************************************************************************/
SamplingSolution: public boolean previous(RunData runData, PhysData physData){
SamplingSolution: return false;
SamplingSolution: } // previous
SamplingSolution:
SamplingSolution: /** Discretize the initial Shape function and initialize the moments
SamplingSolution: Initializes a new set of particles, the number being determined by
SamplingSolution: the numerical scheme parameter.
SamplingSolution: @param function The initial shape to be approximated
SamplingSolution: @see Solution#setScheme
SamplingSolution: ****************************************************************************/
SamplingSolution: public void discretize(ShapeFunction function){
SamplingSolution: int j,k;
SamplingSolution: for (j=0; j<mesh.size(); j++){ //Initial condition
SamplingSolution: f[j] = function.getValue(mesh, j);
SamplingSolution: f0[j] = f[j];
SamplingSolution: g[j] = 0.;
SamplingSolution: }
SamplingSolution:
SamplingSolution: initialMoments[0]=calculateMoments(0); //Conserved quantities
SamplingSolution: initialMoments[1]=calculateMoments(1);
SamplingSolution: initialMoments[2]=calculateMoments(2);
SamplingSolution:
SamplingSolution: for (k=0; k<numberOfRealisations; k++) //Initial position of samplers
SamplingSolution: currentState[k]=function.getSampling(mesh,j);
SamplingSolution:
SamplingSolution: solutionUpToDate = false;
SamplingSolution: } // discretize
SamplingSolution:
SamplingSolution: /** Calculates the deviation from the m:th initial moment.
SamplingSolution: @param m The order of the moment
SamplingSolution: @return The deviation of the m:th moment from the beginning
SamplingSolution: ****************************************************************************/
SamplingSolution: public double momentsDeviation(int m){
SamplingSolution: if(!solutionUpToDate)
SamplingSolution: expectedRealisation();
SamplingSolution: return super.momentsDeviation(m);
SamplingSolution: } // moments
SamplingSolution:
SamplingSolution: /** Calculates the limits of the solution.
SamplingSolution: @return A vector consisting of {min(solution), max(solution)}
SamplingSolution: ****************************************************************************/
SamplingSolution: public double[] limits(){
SamplingSolution: if(!solutionUpToDate)
SamplingSolution: expectedRealisation();
SamplingSolution: return super.limits();
SamplingSolution: } // limits
SamplingSolution:
SamplingSolution: /** Gives the value of the function for an index
SamplingSolution: @param index The index for which to get the value
SamplingSolution: @return The value of the distribution function at a given index
SamplingSolution: ****************************************************************************/
SamplingSolution: public double getValue(int index) {
SamplingSolution: if(!solutionUpToDate)
SamplingSolution: expectedRealisation();
SamplingSolution: return super.getValue(index);
SamplingSolution: } // getValue
SamplingSolution:
SamplingSolution: /** Linear interpolation of the solution. Assumes a uniform mesh.
SamplingSolution: @param arg Argument
SamplingSolution: @return A linear interpolation of the function for a given argument
SamplingSolution: ****************************************************************************/
SamplingSolution: public double getValue(double arg){
SamplingSolution: if(!solutionUpToDate)
SamplingSolution: expectedRealisation();
SamplingSolution: return super.getValue(arg);
SamplingSolution: } // getValue
SamplingSolution:
SamplingSolution: /** Generates a probability distribution and an expectation from the
SamplingSolution: set of sampled values. Assumes a uniform mesh.
SamplingSolution: ****************************************************************************/
SamplingSolution: protected void expectedRealisation(){
SamplingSolution: int j,k;
SamplingSolution: if(pde.equals(jbone.BSCHOLES)){
SamplingSolution: for (j=0; j<mesh.size(); j++){
SamplingSolution: f[j]=0;
SamplingSolution: for (k=0; k<numberOfRealisations; k++) //Arithmetic average estimatate
SamplingSolution: f[j]+= Math.max(currentState[k]-mesh.point(j),0.);
SamplingSolution: f[j]=f[j]/numberOfRealisations;
SamplingSolution: } // for
SamplingSolution: }
SamplingSolution: solutionUpToDate = true;
SamplingSolution:
SamplingSolution: }//expectedRealisation
SamplingSolution:
SamplingSolution: } // SamplingSolution
ShapeBox: /* $Id: ShapeBox.java,v 1.4 2000/05/05 17:32:36 andrej Exp $ */
ShapeBox:
ShapeBox:
ShapeBox: /******************************************************************************
ShapeBox: ShapeBox -- A box function
ShapeBox: @see Solution#discretize
ShapeBox: ******************************************************************************/
ShapeBox: public class ShapeBox extends ShapeFunction{
ShapeBox:
ShapeBox: /** Creates an instance of the class
ShapeBox: @param amplitude The box height
ShapeBox: @param position The box center abscisa
ShapeBox: @param width The box width
ShapeBox: */
ShapeBox: public ShapeBox(double amplitude, double position, double width){
ShapeBox: amplitude_ = amplitude;
ShapeBox: position_ = position;
ShapeBox: width_ = width;
ShapeBox: } // ShapeBox
ShapeBox:
ShapeBox:
ShapeBox: /** Test if the function is complex
ShapeBox: @return true if complex function
ShapeBox: */
ShapeBox: public boolean isComplex(){ return false; };
ShapeBox:
ShapeBox:
ShapeBox: /** Evaluates the shape on a discrete mesh point position.
ShapeBox: @param mesh The mesh coordinates
ShapeBox: @param index The mesh index
ShapeBox: @return The value of the function at the specified mesh location
ShapeBox: */
ShapeBox: public double getValue(Mesh mesh, int index){
ShapeBox: if(mesh.point(index) - position_ > -width_ &&
ShapeBox: mesh.point(index) - position_ < width_)
ShapeBox: return amplitude_;
ShapeBox: else
ShapeBox: return 0;
ShapeBox: } // getValue
ShapeBox:
ShapeBox:
ShapeBox: /** Evaluates complex shape on a discrete mesh point position.
ShapeBox: @param mesh The mesh coordinates
ShapeBox: @param index The mesh index
ShapeBox: @return The value of the function at the specified mesh location
ShapeBox: */
ShapeBox: public Complex getValueC(Mesh mesh, int index){
ShapeBox: Complex z = new Complex(getValue(mesh,index));
ShapeBox: return z;
ShapeBox: } // getValueC
ShapeBox:
ShapeBox:
ShapeBox: /** Returns initial distribution of samplers
ShapeBox: @param mesh The mesh coordinates
ShapeBox: @param index The mesh index
ShapeBox: @return The initial distribution of samplers
ShapeBox: */
ShapeBox: public double getSampling(Mesh mesh, int index){
ShapeBox: double position = position_;
ShapeBox: return position;
ShapeBox: } // getSampling
ShapeBox:
ShapeBox: } // class ShapeBox
ShapeCosinus: /* $Id: ShapeCosinus.java,v 1.4 2000/05/05 17:32:36 andrej Exp $ */
ShapeCosinus:
ShapeCosinus: /******************************************************************************
ShapeCosinus: ShapeCosinus -- A cosinusoidal function
ShapeCosinus: @see Solution#discretize
ShapeCosinus: ******************************************************************************/
ShapeCosinus: public class ShapeCosinus extends ShapeFunction{
ShapeCosinus:
ShapeCosinus: /** Creates an instance of the class
ShapeCosinus: @param amplitude The wave amplitude
ShapeCosinus: @param position The wave phase
ShapeCosinus: @param wavelength The wavelength
ShapeCosinus: */
ShapeCosinus: public ShapeCosinus(double amplitude, double position, double wavelength){
ShapeCosinus: amplitude_ = amplitude;
ShapeCosinus: position_ = position;
ShapeCosinus: wavelength_ = wavelength;
ShapeCosinus: } // ShapeCosinus
ShapeCosinus:
ShapeCosinus:
ShapeCosinus: /** Test if the function is complex
ShapeCosinus: @return true if complex function
ShapeCosinus: */
ShapeCosinus: public boolean isComplex(){ return false; };
ShapeCosinus:
ShapeCosinus:
ShapeCosinus: /** Evaluates the shape on a discrete mesh point position.
ShapeCosinus: @param mesh The mesh coordinates
ShapeCosinus: @param index The mesh index
ShapeCosinus: @return The value of the function at the specified mesh location
ShapeCosinus: */
ShapeCosinus: public double getValue(Mesh mesh, int index){
ShapeCosinus: int n = mesh.size();
ShapeCosinus: double dx = mesh.interval(0);
ShapeCosinus: double a = 2*Math.PI * (double)(index)*dx / wavelength_;
ShapeCosinus: return amplitude_ * (1 + Math.cos(a)) / 1.0;
ShapeCosinus: } // getValue
ShapeCosinus:
ShapeCosinus:
ShapeCosinus: /** Evaluates complex shape on a discrete mesh point position.
ShapeCosinus: @param mesh The mesh coordinates
ShapeCosinus: @param index The mesh index
ShapeCosinus: @return The value of the function at the specified mesh location
ShapeCosinus: */
ShapeCosinus: public Complex getValueC(Mesh mesh, int index){
ShapeCosinus: Complex z = new Complex(getValue(mesh,index));
ShapeCosinus: return z;
ShapeCosinus: } // getValueC
ShapeCosinus:
ShapeCosinus:
ShapeCosinus: /** Returns initial distribution of samplers
ShapeCosinus: @param mesh The mesh coordinates
ShapeCosinus: @param index The mesh index
ShapeCosinus: @return The initial distribution of samplers
ShapeCosinus: */
ShapeCosinus: public double getSampling(Mesh mesh, int index){
ShapeCosinus: double position = position_;
ShapeCosinus: return position;
ShapeCosinus: } // getSampling
ShapeCosinus:
ShapeCosinus: } // class ShapeCosinus
ShapeFunction: /* $Id: ShapeFunction.java,v 1.3 2000/05/05 17:32:36 andrej Exp $ */
ShapeFunction:
ShapeFunction: /** Abstract class for a describing initial shapes
ShapeFunction: @see Solution#discretize
ShapeFunction: */
ShapeFunction: abstract public class ShapeFunction{
ShapeFunction: /** The shape's vertical scale */ protected double amplitude_;
ShapeFunction: /** The shape's horizontal scale */ protected double width_;
ShapeFunction: /** The shape's abscissa */ protected double position_;
ShapeFunction: /** The shape's shorter scale */ protected double wavelength_;
ShapeFunction:
ShapeFunction:
ShapeFunction: /** Test if the function is complex
ShapeFunction: @return true if complex function
ShapeFunction: */
ShapeFunction: abstract public boolean isComplex();
ShapeFunction:
ShapeFunction:
ShapeFunction: /** Evaluates the shape on a discrete mesh point position.
ShapeFunction: @param mesh The mesh coordinates
ShapeFunction: @param index The mesh index
ShapeFunction: @return The value or squared norm of the function at the mesh location
ShapeFunction: */
ShapeFunction: abstract public double getValue(Mesh mesh, int index);
ShapeFunction:
ShapeFunction:
ShapeFunction: /** Evaluates the shape on a discrete mesh point position.
ShapeFunction: @param mesh The mesh coordinates
ShapeFunction: @param index The mesh index
ShapeFunction: @return The value of the function at the specified mesh location
ShapeFunction: */
ShapeFunction: abstract public Complex getValueC(Mesh mesh, int index);
ShapeFunction:
ShapeFunction: /** Returns initial distribution of samplers
ShapeFunction: @param mesh The mesh coordinates
ShapeFunction: @param index The mesh index
ShapeFunction: @return The initial distribution of samplers
ShapeFunction: */
ShapeFunction: abstract public double getSampling(Mesh mesh, int index);
ShapeFunction:
ShapeFunction: } // class ShapeFunction
ShapeFunction:
ShapeGaussian: /* $Id: ShapeGaussian.java,v 1.3 2000/05/05 17:32:36 andrej Exp $ */
ShapeGaussian:
ShapeGaussian: /******************************************************************************
ShapeGaussian: ShapeGaussian -- A Gaussian function
ShapeGaussian: @see Solution#discretize
ShapeGaussian: ******************************************************************************/
ShapeGaussian: public class ShapeGaussian extends ShapeFunction{
ShapeGaussian:
ShapeGaussian: /** Creates an instance of the class
ShapeGaussian: @param amplitude The Gaussian maximum amplitude
ShapeGaussian: @param position The Gaussian abscissa
ShapeGaussian: @param width The Gaussian width
ShapeGaussian: */
ShapeGaussian: public ShapeGaussian(double amplitude, double position, double width){
ShapeGaussian: amplitude_ = amplitude;
ShapeGaussian: position_ = position;
ShapeGaussian: width_ = width;
ShapeGaussian: } // ShapeGaussian
ShapeGaussian:
ShapeGaussian:
ShapeGaussian: /** Test if the function is complex
ShapeGaussian: @return true if complex function
ShapeGaussian: */
ShapeGaussian: public boolean isComplex(){ return false; };
ShapeGaussian:
ShapeGaussian:
ShapeGaussian: /** Evaluates the shape on a discrete mesh point position.
ShapeGaussian: @param mesh The mesh coordinates
ShapeGaussian: @param index The mesh index
ShapeGaussian: @return The value of the function at the specified mesh location
ShapeGaussian: */
ShapeGaussian: public double getValue(Mesh mesh, int index){
ShapeGaussian: double a = (mesh.point(index) - position_) / width_;
ShapeGaussian: return amplitude_ * Math.exp(-a * a);
ShapeGaussian: } // getValue
ShapeGaussian:
ShapeGaussian:
ShapeGaussian: /** Evaluates complex shape on a discrete mesh point position.
ShapeGaussian: @param mesh The mesh coordinates
ShapeGaussian: @param index The mesh index
ShapeGaussian: @return The value of the function at the specified mesh location
ShapeGaussian: */
ShapeGaussian: public Complex getValueC(Mesh mesh, int index){
ShapeGaussian: Complex z = new Complex(getValue(mesh,index));
ShapeGaussian: return z;
ShapeGaussian: } // getValueC
ShapeGaussian:
ShapeGaussian:
ShapeGaussian: /** Returns initial distribution of samplers
ShapeGaussian: @param mesh The mesh coordinates
ShapeGaussian: @param index The mesh index
ShapeGaussian: @return The initial distribution of samplers
ShapeGaussian: */
ShapeGaussian: public double getSampling(Mesh mesh, int index){
ShapeGaussian: double position = position_;
ShapeGaussian: return position;
ShapeGaussian: } // getSampling
ShapeGaussian:
ShapeGaussian: } // class ShapeGaussian
ShapeNumerical: /* $Id: ShapeNumerical.java,v 1.3 2000/05/05 17:32:36 andrej Exp $ */
ShapeNumerical:
ShapeNumerical: /******************************************************************************
ShapeNumerical: ShapeNumerical -- A function (interpolated) from numerical data
ShapeNumerical: @see Solution#discretize
ShapeNumerical: ******************************************************************************/
ShapeNumerical: public class ShapeNumerical extends ShapeFunction{
ShapeNumerical:
ShapeNumerical: /** Numerical function abscissa */ private Mesh myMesh;
ShapeNumerical: /** Numerical function values */ private Solution mySolution;
ShapeNumerical:
ShapeNumerical: /** Creates an instance of the class
ShapeNumerical: @param solution Values defining the shape numerically
ShapeNumerical: */
ShapeNumerical: public ShapeNumerical(Solution solution){
ShapeNumerical: mySolution = solution;
ShapeNumerical: } // ShapeNumerical
ShapeNumerical:
ShapeNumerical:
ShapeNumerical: /** Test if the function is complex
ShapeNumerical: @return true if complex function
ShapeNumerical: */
ShapeNumerical: public boolean isComplex(){ return false; };
ShapeNumerical:
ShapeNumerical:
ShapeNumerical: /** Returns the shape on a discrete mesh point position.
ShapeNumerical: In a future version, this should include a real interpolation,
ShapeNumerical: since the mesh coordinates may not coincide with the data !
ShapeNumerical: @param mesh The mesh coordinates. Now Ignored.
ShapeNumerical: @param index The mesh index
ShapeNumerical: @return The value of the function at the specified mesh location
ShapeNumerical: */
ShapeNumerical: public double getValue(Mesh mesh, int index){
ShapeNumerical: return mySolution.getValue(index);
ShapeNumerical: } // getValue
ShapeNumerical:
ShapeNumerical:
ShapeNumerical: /** Evaluates complex shape on a discrete mesh point position.
ShapeNumerical: @param mesh The mesh coordinates
ShapeNumerical: @param index The mesh index
ShapeNumerical: @return The value of the function at the specified mesh location
ShapeNumerical: */
ShapeNumerical: public Complex getValueC(Mesh mesh, int index){
ShapeNumerical: Complex z = new Complex(getValue(mesh,index));
ShapeNumerical: return z;
ShapeNumerical: } // getValueC
ShapeNumerical:
ShapeNumerical:
ShapeNumerical: /** Returns initial distribution of samplers
ShapeNumerical: @param mesh The mesh coordinates
ShapeNumerical: @param index The mesh index
ShapeNumerical: @return The initial distribution of samplers
ShapeNumerical: */
ShapeNumerical: public double getSampling(Mesh mesh, int index){
ShapeNumerical: double position = position_;
ShapeNumerical: return position;
ShapeNumerical: } // getSampling
ShapeNumerical:
ShapeNumerical: } // class ShapeNumerical
ShapeOption: /* $Id: ShapeOption.java,v 1.1 2000/05/05 17:32:36 andrej Exp $ */
ShapeOption:
ShapeOption: /******************************************************************************
ShapeOption: ShapeOption -- Payoff function of an option at expiry
ShapeOption: @see Solution#discretize
ShapeOption: ******************************************************************************/
ShapeOption: public class ShapeOption extends ShapeFunction{
ShapeOption:
ShapeOption: /** Creates an instance of the class
ShapeOption: @param position The exercise price (strike)
ShapeOption: @param amplitude Not used so far
ShapeOption: @param width Not used so far
ShapeOption: */
ShapeOption: public ShapeOption(double amplitude, double position, double width){
ShapeOption: position_ = position;
ShapeOption: amplitude_ = amplitude;
ShapeOption: width_ = width;
ShapeOption: } // ShapeOption
ShapeOption:
ShapeOption:
ShapeOption: /** Test if the function is complex
ShapeOption: @return true if complex function
ShapeOption: */
ShapeOption: public boolean isComplex(){ return false; };
ShapeOption:
ShapeOption:
ShapeOption: /** Evaluates the shape on a discrete mesh point position.
ShapeOption: @param mesh The mesh coordinates
ShapeOption: @param index The mesh index
ShapeOption: @return The value of the function at the specified mesh location
ShapeOption: */
ShapeOption: public double getValue(Mesh mesh, int index){
ShapeOption: double value = Math.max(position_ - mesh.point(index),0.);
ShapeOption: return value;
ShapeOption: } // getValue
ShapeOption:
ShapeOption:
ShapeOption: /** Evaluates complex shape on a discrete mesh point position.
ShapeOption: @param mesh The mesh coordinates
ShapeOption: @param index The mesh index
ShapeOption: @return The value of the function at the specified mesh location
ShapeOption: */
ShapeOption: public Complex getValueC(Mesh mesh, int index){
ShapeOption: Complex z = new Complex(getValue(mesh,index));
ShapeOption: return z;
ShapeOption: } // getValueC
ShapeOption:
ShapeOption:
ShapeOption: /** Returns initial distribution of samplers
ShapeOption: @param mesh The mesh coordinates
ShapeOption: @param index The mesh index
ShapeOption: @return The initial distribution of samplers
ShapeOption: */
ShapeOption: public double getSampling(Mesh mesh, int index){
ShapeOption: double position = position_;
ShapeOption: return position;
ShapeOption: } // getSampling
ShapeOption:
ShapeOption: } // class ShapeOption
ShapeSoliton: /* $Id: ShapeSoliton.java,v 1.5 2000/05/05 17:32:36 andrej Exp $ */
ShapeSoliton:
ShapeSoliton: /******************************************************************************
ShapeSoliton: ShapeSoliton -- A soliton function
ShapeSoliton: @see Solution#discretize
ShapeSoliton: ******************************************************************************/
ShapeSoliton: public class ShapeSoliton extends ShapeFunction{
ShapeSoliton:
ShapeSoliton: /** Creates an instance of the class
ShapeSoliton: @param amplitude The soliton amplitude
ShapeSoliton: @param position The soliton position
ShapeSoliton: @param width Is not appropriate for this shape
ShapeSoliton: */
ShapeSoliton: public ShapeSoliton(double amplitude, double position, double width){
ShapeSoliton: amplitude_ = amplitude;
ShapeSoliton: position_ = position;
ShapeSoliton: width_ = width;
ShapeSoliton: } // ShapeSoliton
ShapeSoliton:
ShapeSoliton:
ShapeSoliton: /** Test if the function is complex
ShapeSoliton: @return true if complex function
ShapeSoliton: */
ShapeSoliton: public boolean isComplex(){ return false; };
ShapeSoliton:
ShapeSoliton:
ShapeSoliton: /** Evaluates the shape on a discrete mesh point position.
ShapeSoliton: @param mesh The mesh coordinates
ShapeSoliton: @param index The mesh index
ShapeSoliton: @return The value of the function at the specified mesh location
ShapeSoliton: */
ShapeSoliton: public double getValue(Mesh mesh, int index){
ShapeSoliton: int n = mesh.size();
ShapeSoliton: double dx = mesh.interval(0);
ShapeSoliton: double len = (double)n*dx;
ShapeSoliton: double sca = 2./3.*amplitude_;
ShapeSoliton: double a = (mesh.point(index) -
ShapeSoliton: (len-position_))*Math.sqrt(sca/2);
ShapeSoliton: double result = 3*sca/Math.pow( 0.5*(Math.exp(a)+Math.exp(-a)), 2);
ShapeSoliton: sca=amplitude_ / 3.0;
ShapeSoliton: a=(mesh.point(index)-position_) *
ShapeSoliton: Math.sqrt(sca/2);
ShapeSoliton: result += 3*sca/Math.pow( 0.5*(Math.exp(a)+Math.exp(-a)), 2);
ShapeSoliton: return result;
ShapeSoliton: } // getValue
ShapeSoliton:
ShapeSoliton:
ShapeSoliton: /** Evaluates complex shape on a discrete mesh point position.
ShapeSoliton: @param mesh The mesh coordinates
ShapeSoliton: @param index The mesh index
ShapeSoliton: @return The value of the function at the specified mesh location
ShapeSoliton: */
ShapeSoliton: public Complex getValueC(Mesh mesh, int index){
ShapeSoliton: Complex z = new Complex(getValue(mesh,index));
ShapeSoliton: return z;
ShapeSoliton: } // getValueC
ShapeSoliton:
ShapeSoliton:
ShapeSoliton: /** Returns initial distribution of samplers
ShapeSoliton: @param mesh The mesh coordinates
ShapeSoliton: @param index The mesh index
ShapeSoliton: @return The initial distribution of samplers
ShapeSoliton: */
ShapeSoliton: public double getSampling(Mesh mesh, int index){
ShapeSoliton: double position = position_;
ShapeSoliton: return position;
ShapeSoliton: } // getSampling
ShapeSoliton:
ShapeSoliton: } // class ShapeSoliton
ShapeWavePacket: /* $Id: ShapeWavePacket.java,v 1.2 2000/05/05 17:32:37 andrej Exp $ */
ShapeWavePacket:
ShapeWavePacket: /******************************************************************************
ShapeWavePacket: ShapeWavePacket -- A Wave Packet function
ShapeWavePacket: @see Solution#discretize
ShapeWavePacket: ******************************************************************************/
ShapeWavePacket: public class ShapeWavePacket extends ShapeFunction{
ShapeWavePacket:
ShapeWavePacket: /** Creates an instance of the class
ShapeWavePacket: @param amplitude The wave packet amplitude
ShapeWavePacket: @param width The localization
ShapeWavePacket: @param position The absissa
ShapeWavePacket: */
ShapeWavePacket: public ShapeWavePacket(double amplitude, double position, double width,
ShapeWavePacket: double wavelength){
ShapeWavePacket: amplitude_ = amplitude;
ShapeWavePacket: position_ = position;
ShapeWavePacket: width_ = width;
ShapeWavePacket: wavelength_ = wavelength;
ShapeWavePacket: } // ShapeWavePacket
ShapeWavePacket:
ShapeWavePacket:
ShapeWavePacket: /** Test if the function is complex
ShapeWavePacket: @return true if complex function
ShapeWavePacket: */
ShapeWavePacket: public boolean isComplex(){ return true; };
ShapeWavePacket:
ShapeWavePacket:
ShapeWavePacket: /** Evaluates the shape on a discrete mesh point position.
ShapeWavePacket: @param mesh The mesh coordinates
ShapeWavePacket: @param index The mesh index
ShapeWavePacket: @return The value of the function at the specified mesh location
ShapeWavePacket: */
ShapeWavePacket: public Complex getValueC(Mesh mesh, int index){
ShapeWavePacket: double norm = Math.sqrt(Math.PI/width_);
ShapeWavePacket: double x0 = (mesh.point(index) - position_);
ShapeWavePacket: Complex osc = new Complex(0.,2.*Math.PI/wavelength_*x0);
ShapeWavePacket: Complex psi = new Complex(osc.exp().scale(
ShapeWavePacket: amplitude_*Math.exp(-x0*x0/(4.*width_))));
ShapeWavePacket: return psi;
ShapeWavePacket:
ShapeWavePacket: } // getValueC
ShapeWavePacket:
ShapeWavePacket: /** Evaluates the square of the norm
ShapeWavePacket: @param mesh The mesh coordinates
ShapeWavePacket: @param index The mesh index
ShapeWavePacket: @return The value of the function at the specified mesh location
ShapeWavePacket: */
ShapeWavePacket: public double getValue(Mesh mesh, int index){
ShapeWavePacket: Complex psi = getValueC(mesh,index);
ShapeWavePacket: double z = psi.norm();
ShapeWavePacket: return z;
ShapeWavePacket:
ShapeWavePacket: } // getValue
ShapeWavePacket:
ShapeWavePacket:
ShapeWavePacket: /** Returns initial distribution of samplers
ShapeWavePacket: @param mesh The mesh coordinates
ShapeWavePacket: @param index The mesh index
ShapeWavePacket: @return The initial distribution of samplers
ShapeWavePacket: */
ShapeWavePacket: public double getSampling(Mesh mesh, int index){
ShapeWavePacket: double position = position_;
ShapeWavePacket: return position;
ShapeWavePacket: } // getSampling
ShapeWavePacket:
ShapeWavePacket: } // class ShapeWavePacket
Solution: /* $Id: Solution.java,v 1.38 2001/03/24 20:23:56 andrej Exp $ */
Solution:
Solution: import java.awt.*;
Solution: import java.text.*;
Solution:
Solution: /******************************************************************************
Solution: * Solution -- Is an abstract class containing all the solutions and solvers.
Solution: * @see FluidSolution
Solution: * @see ParticleSolution
Solution: ******************************************************************************/
Solution: abstract public class Solution{
Solution:
Solution: /** The PDE to be solved @see #setPde */
Solution: protected String pde;
Solution: /** Numerical method name @see #setMethod */
Solution: protected String method;
Solution: /** Numerical scheme name @see #setScheme */
Solution: protected String scheme;
Solution: /** Initial condition name @see #setIC */
Solution: protected String ic;
Solution: /** Absolute time @see #setTime*/
Solution: protected double time;
Solution:
Solution: /** Function time level 0 (IC) @see #f */ protected double[] f0;
Solution:
Solution: /** Function time level n-1 */ protected double[] fm;
Solution: /** Function time level n */ protected double[] f;
Solution: /** Function time level n+1 */ protected double[] fp;
Solution: /** Derivative time level n-1 @see #fm */ protected double[] dfm;
Solution: /** Derivative time level n @see #f */ protected double[] df;
Solution: /** Derivative time level n+1 @see #fp */ protected double[] dfp;
Solution:
Solution: /** Extra Function time level n-1 */ protected double[] gm;
Solution: /** Extra Function time level n */ protected double[] g;
Solution: /** Extra Function time level n+1 */ protected double[] gp;
Solution: /** Extra Derivative time level n-1 @see #gm */protected double[] dgm;
Solution: /** Extra Derivative time level n @see #g */ protected double[] dg;
Solution: /** Extra Derivative time level n+1 @see #gp */protected double[] dgp;
Solution:
Solution: /** Complex Function time level n-1 */ protected Complex[] hm;
Solution: /** Complex Function time level n */ protected Complex[] h;
Solution: /** Complex Function time level n+1 */ protected Complex[] hp;
Solution: /** Derivative time level n-1 @see #hm */ protected Complex[] dhm;
Solution: /** Derivative time level n @see #h */ protected Complex[] dh;
Solution: /** Derivative time level n+1 @see #hp */ protected Complex[] dhp;
Solution:
Solution: /** Complex/spectrum time level n-1 */ protected Complex[] sm;
Solution: /** Complex/spectrum time level n */ protected Complex[] s;
Solution: /** Complex/spectrum time level n+1 */ protected Complex[] sp;
Solution:
Solution: /** The underlaying mesh of the solution */ protected Mesh mesh;
Solution: /** The mesh points @see #mesh */ protected double[] x;
Solution: /** The mesh intervals @see #mesh */ protected double[] dx;
Solution: /** Store initial moments */ protected double[]
Solution: initialMoments = {0.,0.,0.};
Solution:
Solution: /** Horizontal scale of plot area */ protected int xSize;
Solution: /** Vertical scale of plot area */ protected int ySize;
Solution: /** Horizontal offset of plot area */ protected int xOffset;
Solution: /** Vertical offset of plot area */ protected int yOffset;
Solution:
Solution: /** Creates an instance of a Solution object with all the attributes.
Solution: @param runData The run time parameters
Solution: ****************************************************************************/
Solution: public Solution(RunData runData) {
Solution: Mesh mesh = runData.getMesh();
Solution: this.mesh = mesh;
Solution: int n = mesh.size();
Solution: x = mesh.points(); dx = mesh.intervals(); f0= new double[n];
Solution: fm = new double[n]; f = new double[n]; fp = new double[n];
Solution: dfm= new double[n]; df = new double[n]; dfp = new double[n];
Solution: gm = new double[n]; g = new double[n]; gp = new double[n];
Solution: dgm= new double[n]; dg = new double[n]; dgp = new double[n];
Solution: hm = new Complex[n]; h = new Complex[n]; hp = new Complex[n];
Solution: dhm= new Complex[n]; dh = new Complex[n]; dhp = new Complex[n];
Solution: sm = new Complex[n]; s = new Complex[n]; sp = new Complex[n];
Solution: lastStep = 0;
Solution: lastTimeStep = 0;
Solution: } // Solution
Solution:
Solution:
Solution: /** Internal function to calculate the initial moments of f[] or their
Solution: deviation from the beginning of the evolution.
Solution: @param m The order of the moment
Solution: @return The value or deviation of the m:th moment of f[]
Solution: @see #momentsDeviation
Solution: ****************************************************************************/
Solution: protected double calculateMoments(int m){
Solution: int n=f.length-1;
Solution: double moment = 0.0;
Solution: if (m==0) {
Solution: if (pde.equals(jbone.EXERCISE) &&
Solution: scheme.equals(jbone.EXC3_02)) { //Cylindrical Jacobian
Solution: for (int i=0; i<n; i++) {
Solution: moment= moment+0.5*dx[i]*(f[i+1]+f[i])*0.5*(x[i]+x[i+1]);
Solution: }
Solution: } else { //Cartesian Jacobian
Solution: for (int i=0; i<n; i++) {
Solution: moment= moment+0.5*dx[i]*(f[i+1]+f[i]);
Solution: }
Solution: moment = moment +0.5*dx[n]*(f[0]+f[n]);
Solution: }
Solution: } else { moment=0.0; }
Solution: return moment; //Initial value
Solution: } // calculateMoments
Solution:
Solution: /** Calculates the deviation from the m:th initial moment.
Solution: @param m The order of the moment
Solution: @return The deviation of the m:th moment from the beginning
Solution: ****************************************************************************/
Solution: public double momentsDeviation(int m){
Solution: return (calculateMoments(m) - initialMoments[m])
Solution: / initialMoments[m];
Solution: } // moments
Solution:
Solution: /** Calculates the solution boundaries.
Solution: @return A vector with {min(solution), max(solution)}
Solution: ****************************************************************************/
Solution: protected double[] limits() {
Solution: double[] lim = {f[0], f[0]};
Solution: for (int i=1;i<f.length;i++) {
Solution: if (f[i]<lim[0]) {lim[0]=f[i];};
Solution: if (f[i]>lim[1]) {lim[1]=f[i];};
Solution: if (g[i]!=0. && g[i]<lim[0]) {lim[0]=g[i];};
Solution: if (g[i]!=0. && g[i]>lim[1]) {lim[1]=g[i];};
Solution: } // for
Solution: return lim;
Solution: } // limits
Solution:
Solution:
Solution: /** Set the PDE as a string.
Solution: @param pde The name of the equation
Solution: @see jbone#ADVECTION
Solution: ****************************************************************************/
Solution: public void setPde(String pde){
Solution: this.pde = new String(pde);
Solution: } // setPde
Solution:
Solution: /** Set the numerical method as a string.
Solution: @param method The name of the numerical method
Solution: @return True
Solution: ****************************************************************************/
Solution: public void setMethod(String method){
Solution: this.method = new String(method);
Solution: } // setMethod
Solution:
Solution: /** Set the numerical scheme as a string.
Solution: @param scheme The name of the numerical scheme
Solution: @return True
Solution: ****************************************************************************/
Solution: public void setScheme(String scheme){
Solution: this.scheme = new String(scheme);
Solution: } // setScheme
Solution:
Solution: /** Set the initial condision as a string.
Solution: @param ic The name of the initial condition
Solution: @return True
Solution: ****************************************************************************/
Solution: public void setIC(String ic){
Solution: this.ic = new String(ic);
Solution: } // setIC
Solution:
Solution:
Solution: /** Set the current physical time
Solution: @param time Current time
Solution: @return True
Solution: ****************************************************************************/
Solution: public boolean setTime(double time){
Solution: this.time = time;
Solution: return true;
Solution: } // setTime
Solution:
Solution:
Solution: /** Increment the current physical time
Solution: @param time Current time increment
Solution: @return True
Solution: ****************************************************************************/
Solution: public boolean incTime(double time){
Solution: this.time += time;
Solution: return true;
Solution: } // incTime
Solution:
Solution:
Solution: /** Get the current physical time
Solution: @return current time
Solution: ****************************************************************************/
Solution: public double getTime(){
Solution: return this.time;
Solution: } // getTime
Solution:
Solution:
Solution: /** Discretize the initial shape and initialize the moments
Solution: @param function The initial shape to be approximated
Solution: ****************************************************************************/
Solution: abstract public void discretize(ShapeFunction function);
Solution:
Solution: /** The current step in the time discretization */
Solution: private int lastStep;
Solution: /** The last time step used */
Solution: private double lastTimeStep;
Solution:
Solution: /** Updates the information for the headers
Solution: @param runData List of run parameters
Solution: @param step The current step in the simulation
Solution: @see #plot
Solution: ****************************************************************************/
Solution: public void updateHeaders(RunData runData, int step){
Solution: lastStep = step;
Solution: lastTimeStep = runData.getParamValue(runData.timeStepNm);
Solution: } // updateHeaders
Solution:
Solution: /** Advance the solution forward one step in time.
Solution: The equation is set by the internal variable equation_ and
Solution: the numerical scheme by the internal variable scheme_
Solution: @param runData List of run parameters
Solution: @param physData Physical parameters (e.g. potential)
Solution: @see RunData
Solution: @see PhysData
Solution: @return True if a scheme is implemented for that equation
Solution: */
Solution: abstract public boolean next(RunData runData, PhysData physData);
Solution:
Solution: /** Take the solution backward one step to initialize schemes
Solution: with 3 time levels.
Solution: The equation is set by the internal variable equation_ and
Solution: the numerical scheme by the internal variable scheme_
Solution: @param runData List of run parameters
Solution: @see RunData
Solution: @return True if an initialisation scheme is implemented for that equation
Solution: */
Solution: abstract public boolean previous(RunData runData, PhysData physData);
Solution:
Solution: /** Gives the value of the function for an index
Solution: @param index The index for which to get the value
Solution: @return The value of the function at a given index
Solution: ****************************************************************************/
Solution: public double getValue(int index) {
Solution: return f[index];
Solution: } // getValue
Solution:
Solution: /** Linear interpolation of the solution. Assumes a uniform mesh.
Solution: @param arg Argument
Solution: @return A linear interpolation of the function for a given argument
Solution: ****************************************************************************/
Solution: public double getValue(double arg){
Solution: double[] xLim = mesh.limits();
Solution: int index=(int)Math.floor((arg-xLim[0])* // Lower mesh cell index
Solution: (mesh.size()-1)/(xLim[1]-xLim[0]));
Solution: if(index == mesh.size()-1) index--;
Solution: double x_low = mesh.point(index); //Linear interpolation
Solution: double x_high = mesh.point(index+1);
Solution: double a = (arg-x_low)/(x_high-x_low);
Solution: return f[index]*(1-a) + f[index+1]*a;
Solution: } // getValue
Solution:
Solution: /** Tells whether the solution implements a option
Solution: @param option The the option to implement
Solution: @return True if the option is implemented
Solution: @see jbone#pdeNames
Solution: @see jbone#schemeNames
Solution: */
Solution: abstract public boolean hasOption(String option);
Solution:
Solution: /** Set the corners of the plot area
Solution: ****************************************************************************/
Solution: public void rescale(){
Solution: double[] lim = limits();
Solution: // Extra space needed for the periodic boundary condition line
Solution: x_0 = x[0] - dx[0] / 2.0;
Solution: x_1 = x[x.length - 1] + dx[x.length - 1] / 2.0;
Solution: y_0 = lim[0];
Solution: y_1 = lim[1];
Solution: if (y_0 == y_1) { y_0=- 1.; y_1=+ 1.; }
Solution: } // rescale
Solution:
Solution: /** x of lower left corner of plot area */
Solution: protected double x_0;
Solution: /** x of upper right corner of plot area */
Solution: protected double x_1;
Solution: /** y of lower left corner of plot area */
Solution: protected double y_0;
Solution: /** y of upper right corner of plot area */
Solution: protected double y_1;
Solution:
Solution: /** Returns window size
Solution: @see #plot
Solution: ****************************************************************************/
Solution: public int[] getWinSize(){
Solution: int[] ret = {xSize,ySize,xOffset,yOffset};
Solution: return ret;
Solution: } //getWinSize
Solution:
Solution: /** Get physical coordinates from mouse coordinates
Solution: @param x horizontal mouse coordinate
Solution: @param y vertical mouse coordinate
Solution: @see #plot
Solution: ****************************************************************************/
Solution: public double[] measure(int x, int y){
Solution: double dx = (xSize-2*xOffset)/(x_1-x_0);
Solution: double dy = (ySize-2*yOffset)/(y_1-y_0);
Solution: double X = x_0+((double)( x-xOffset))/dx;
Solution: double Y = y_0+((double)(-y+ySize+1+5))/dy;
Solution: double[] coord = {X,Y};
Solution: return coord;
Solution: } // measure
Solution:
Solution: /** Plots the solution
Solution: @param plotArea The plot area
Solution: @param offScrImage The off screen image to draw on
Solution: @param headers Whether to draw headers
Solution: ****************************************************************************/
Solution: public void plot(Canvas plotArea, Image offScrImage, boolean headers){
Solution: } // plot
Solution:
Solution: /** Print the solution in ASCII to java console
Solution: ****************************************************************************/
Solution: public void output(int step){
Solution: System.out.println(" ");
Solution: System.out.println("Step="+step);
Solution: for (int i=0; i<mesh.size(); i++) {
Solution: System.out.println("x["+i+"]="+mesh.point(i)+
Solution: " f0["+i+"]="+f0[i]+
Solution: " f["+i+"]="+f[i]);
Solution: }
Solution: } // print
Solution:
Solution: } // class Solution
|
