/* sim3.java contains several CSimulation subclasses */ /* Part of the www.MyPhysicsLab.com physics simulation applet. Copyright (c) 2001 Erik Neumann This program is free software; you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation; either version 2 of the License, or (at your option) any later version. This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received a copy of the GNU General Public License along with this program; if not, write to the Free Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA Contact Erik Neumann at erikn@MyPhysicsLab.com or 610 N. 65th St. Seattle WA 98103 */ import java.applet.*; import java.awt.*; /* o String Gershenfeld, Nature of Mathematical Modeling, gives a way of discretizing the PDE for a string. I'm imagining a string fixed at each end, and you can choose initial conditions somehow (eg. drag at a single point) and then run the simulation. Maybe have an oscillating driver, like a cello bow, that applies a force at a certain point repetitively. # find equations (ie. need tension & mass factors...) # new element that draws string corresponding to the vector. # create the u() vector (maybe 1000 long) # have element sample it at around 100 places for drawing # Test by evolving according to some simple formula. eg. u(x,t) = sin(2*pi*x/length)*sin(t) # 4 large vectors: old slope, old velocity, new slope, new velocity (note: can save time copying over by having a flag that says which is the more recent vector). # can just set initial conditions by code at first The velocity vector can be all zero. Fill in the slope vector by deciding on a function, and calc its first derivative. # solver applies eqns to calculate new from old. # integration of slope vector fills in the u vector. This can be done in doModifyObjects. # stability condition is: abs(v) delta(t) / delta(x) < 1 Let's check what is the delta(t) that we are actually getting in this and other simulations. If its too big, then we need to go to non-real time, or find ways to speed things up! Note that we are very limited by delta(t)... can't choose delta(x) to be very small without also decreasing delta(t). We can artifically slow down time by having a scaling factor on time. Make this a variable that is controllable, and then we can see effect when its non-stable. # Show the stability as a (non-modifiable) number in control window. (or in message window as a quick cludge). # note: had to fix boundary conditions (had wrong assumptions, they are not both zero, only the velocity is zero). # The problem with (nearly) square waves: Seems to get 2 different solutions that oscillate. This is an interesting artifact of the numerical method. Really there are two independent simulations going on, one involves even numbered slots (indices) in the vectors and the other involving odd numbered slots. Let r[t,x] be the r (slope) vector at time t, position x. and similarly for the s (speed) vector. Algorithm is then r[t+1,j] = (1/2)(r[t,j+1] + r[t,j-1]) + k(s[t,j+1] - s[t,j-1]) and similarly for s, where k is a constant. Therefore, the even numbered slots depend only on the odd numbered slots in the last time period, and vice versa. Therefore, if the initial conditions are such that r[0,j] is very different from r[0,j+1], then you wind up with two different simulations for the odd and even slots. In the square wave case, I had all of r & s zero except for two points in r: r[2*n/5] and r[3*n/5]. Therefore, "half" of the simulation saw entirely zero vectors and the other half saw those two points. The reason that the resulting string oscillated is that the integration algorithm adds even & odd spots at different strengths (2:1) and since the non-zero simulation switches between even and odd slots every iteration, this caused the oscillation. Is there perhaps a way to connect the two separate even/odd simulation and keep things together? # Try staggered leapfrog method This method didn't work either.... # What works: Algorithm 12.4, "Wave Equation Finite-Difference" from Numerical Analysis, 6th Ed. by Burden & Faires. This definitely works where the others I tried (from Numerical Recipes) were failing. Actually the Lax method sort of seemed to work, but had a lot of numerical dispersion (smoothing over time). # manual forcing curve element should take up less than full screen (maybe draw a boundary around it?) set a block at one end and allow it to be moved up & down set the endpoint of the string to be position of block. GOT IT! It works! # select starting waveform from a menu Need to add a new menu that is different for each simulation. Is it a single menu which changes items, or different menus? # compare to ANALYTIC solution Plot both, with different colors. o Reset position of block when restarting (ie. option menu select). o Change delta(t) and delta(x) as parameters? Need to show stability somewhere.... o specify amplitude of first 10 harmonics, and their phases? o Reset button o repetitive FORCING: add parameters for amplitude & frequency of forcing, show how this leads to "ringing" when frequency is same as natural frequency of the string. Might need damping to get this. o PULSE Provide a single (half-cycle) of forcing, to create a smooth pulse. Can we solve this analytically? This one seems to build up trailing artifacts over time... why? o Show FREQUENCY ANALYSIS of the string Possible to show frequency breakdown of the string? Maybe could see loss of high-frequency components over time? Is this what causes the build-up of wobbles in a wave over time (and loss of shape of a square wave). o Try using smaller time deltas (this will slow down the simulation) and see if the artifacts are reduced. o Ability to draw the starting waveform, or o Need to understand how the Burden & Faires method is derived. Work it out in Mathematica... then can maybe apply similar results to adding forcing or damping or gravity, or a weight at certain points, etc. o add in gravity or damping or forcing, see Tung p. 29 To do this, need to re-derive the recursion relationships from the new PDE. o Numerical considerations are mainly about how to represent the derivatives numerically in the recursion, and stability. o SPEED how to make the whole thing go faster? Need to increase tension (decrease density), but must keep stability below 1. Therefore, have to decrease delta(t). But we are only getting about 20 frames per second. Where is all the time going? Check in simpler simulations, what is maximum achievable? Try profiling. Use full-screen mode to get more time. */ class CCurve extends CElement { protected static final int DRAW_POINTS = 21; public double[] m_Data = null; // curve is defined by this data vector public double m_Height; public double m_Width; public CCurve (double X1, double Y1, double width, double height) { super(X1, Y1, X1+width, Y1+height); m_Height = height; m_Width = width; } public void draw (Graphics g, CoordMap map) { // Draw the curve defined by the data vector. // Strategy: // There are m_N points in the data vector // We want to sample the data vector at a lesser number (= DRAW_POINTS) // of points. To do so, we use floating point calculations to pick // the x-values and then determine the nearest index to those points // in the data vector. // WARNING: the following assumes DRAW_POINTS <= m_N ???? double x,y,x0,y0; int i,j; if (m_Data == null) return; x0 = 0; y0 = m_Data[0]; int scrx = map.simToScreenX(m_X1+x0); int scry = map.simToScreenY(y0); int oldx, oldy; g.setColor(m_Color); for(i=1;i= m_Data.length) j = m_Data.length-1; //System.out.println("draw point "+j); y = m_Data[j]; oldx = scrx; oldy = scry; scrx = map.simToScreenX(m_X1+x); scry = map.simToScreenY(y); g.drawLine(oldx, oldy, scrx, scry); x0 = x; y0 = y; } } } class CString extends CSimulation { protected static final int FLAT = 0; protected static final int TRIANGLE = 1; protected static final int QUARTER_TRI = 2; protected static final int SQUARE_PULSE = 3; protected static final int SINE_PULSE = 4; protected static final int HALF_SINE_PULSE = 5; protected static final int FANCY_SINE = 6; CCurve m_Curve; CCurve m_Curve2; // analytic version CMass m_Block; protected double m_Length; // length of string; protected double m_Tension; // tension of string protected double m_Density; // density of string per unit length protected double gravity; int m_shape; // starting wave shape private static final int STRING_POINTS = 21; private double[] w1; // data array private double[] w2; // data array private double[] w3; // data array private double[] w4; // analytic data array int w_idx; // tells which array (1 or 2 or 3) is most recent double delta_x; // delta x double delta_t; // delta time (fixed) double total_t; // cumulative time; private static final int AVG_LEN = 10; double[] times; // for checking on delta(t) = avg time between updates double[] stab; // for averaging stability value int time_idx; // index into times[] array double last_time; // last time we printed delta(t) /* set initial conditions for string. input is width of the string */ void Initialize_Shape() { w_idx = 2; w1 = new double[STRING_POINTS]; w2 = new double[STRING_POINTS]; w3 = new double[STRING_POINTS]; w4 = new double[STRING_POINTS]; double r = (delta_t*delta_t*m_Tension/m_Density)/(delta_x*delta_x); w1[0] = w1[STRING_POINTS-1] = 0; w2[0] = w2[STRING_POINTS-1] = 0; int i; for (i=1;i w) return 0; else return w*((x < w/2) ? x/w : 1 - x/w); case SQUARE_PULSE: w = m_Length/4; return (x < w) ? w : 0; case SINE_PULSE: w = m_Length/3; if (x>w) return 0; else return Math.sin(2*Math.PI*x/w); case HALF_SINE_PULSE: w = m_Length/3; if (x>w) return 0; else return Math.sin(Math.PI*x/w); case FANCY_SINE: return (Math.sin(2*Math.PI*x/m_Length) + Math.sin(4*Math.PI*x/m_Length) + Math.sin(6*Math.PI*x/m_Length))/3; } } public CString(Applet app, int shape) { super(app); m_shape = shape; numVars = 0; var_names = new String[] {}; // empty string //FillTestData(pi/2); // initialize stuff for figuring average time delta times = new double[AVG_LEN]; // automatically initialized to zero stab = new double[AVG_LEN]; last_time = 0; time_idx = 0; map = new CoordMap(CoordMap.INCREASE_UP, 0, 15, -2, 2, CoordMap.ALIGN_MIDDLE, CoordMap.ALIGN_MIDDLE); // determine window size in simulation coords, // to make the curve fill the window m_Length = 3*map.sim_width/4; m_Tension = 30; m_Density = 1; gravity = 9; delta_x = m_Length/(STRING_POINTS-1); // delta x // delta_t = 0.055; // delta time is fixed!!! delta_t = 0.03; total_t = 0; Initialize_Shape(); double b = 0.7; // width of block double h = map.sim_height/2; m_Curve = new CCurve(map.sim_x1+2*b, -h, m_Length, 2*h); //X1, Y1, width, height add(m_Curve); // second curve is to show analytic version m_Curve2 = new CCurve(map.sim_x1+2*b, -h, m_Length, 2*h); m_Curve2.m_Data = w4; m_Curve2.m_Color = Color.red; // This mass is for user to drag up & down for manual forcing m_Block = new CMass(m_Curve.m_X1 - b, -b/2, b, b, CElement.MODE_RECT); add(m_Block); doModifyObjects(); params = new String[] {"gravity","tension","density"}; } public double get(int param) { switch(param) { case 0: return gravity; case 1: return m_Tension; case 2: return m_Density; default: return 0; } } public void set(int param, double value) { switch(param) { case 0: gravity = value; break; case 1: m_Tension = value; break; case 2: m_Density = value; break; } } public void doModifyObjects() { // set the curve to point at the current vector double[] w; switch (w_idx) { default: case 1: w = w1; break; case 2: w = w2; break; case 3: w = w3; break; } m_Curve.m_Data = w; } public void startDrag(CElement e) { m_Animating = true; } public void constrainedSet(CElement e, double x, double y) { if (e==m_Block) { e.setY1(y); } } // don't use the Runge Kutta diff eq solver // instead calculate the variables directly from the time. public void solve(double t, double h) { int i; total_t += delta_t; CalcAnalyticData(total_t); //if (t>-99999) // return; // figure out which vector to use for latest data double[] w_new; double[] w; double[] w_old; switch (w_idx) { default: case 1: // w1 is most recent data, then 3, 2 is oldest w = w1; w_old = w3; w_new = w2; w_idx = 2; break; case 2: // w2 is most recent data, then 1, 3 is oldest w = w2; w_old = w1; w_new = w3; w_idx = 3; break; case 3: // w3 is most recent data, then 2, 1 is oldest w = w3; w_old = w2; w_new = w1; w_idx = 1; break; } double r = (delta_t*delta_t*m_Tension/m_Density)/(delta_x*delta_x); // here is the heart of the PDE solver for (i=1; i= AVG_LEN) time_idx = 0; times[time_idx] = h; stab[time_idx] = Math.sqrt(r); if (t - last_time > 2) { last_time = t; double d=0; double e=0; for (i=0; i