Solving Poisson Equation Numerically

Potential between Plates 0, 0, and 5, 5

In the absence of charge, the Poisson equation reduces to the Laplace equation .
The Laplace equation has the nice property that the solution is the average of the surrounding points. This is applicable for the case of one, two, and three dimensions. The reason for this is that the average of the potential on a sphere is given by where q is the charge inside the sphere. Since the region we are considering does not have any volume charge density, this is zero. Thus, the average field is just the field at the center of the sphere. The one dimensional case is not very interesting: it simply states that the function that satisfies the Laplace equation is linear, of the form ax + b. The two dimensional case is more interesting. Usually, the voltage at the boundary is specified, and this value can be used to calculate the value of the voltage everywhere. The program starts with the known values of voltage, and iteratively builds the voltage at the interior points. The averaging scheme used minimizes errors. Of course, I could simply have taken the values above, below, left and right, and that would’ve given the answer too (eventually), but this can be shown mathematically to be the one with least error. Note that this can also be used in the case where the charge density inside is non zero, but now the equation will have to be modified, as it will not be a simple average.

The problem that is considered here can be solved using separation of variables, so that an analytical solution exists. However, it can be used in cases where such methods are not that easily applied.

#include

int main(int argc, char** argv){
float pot[50][50];
float fc, fs;
int i, j, k;

for (i = 0; i < 50; i++) for (j = 0; j < 50; j++) pot[i][j] = 0; for (i = 0;i < 50; i++) { pot[0][i] = 5; pot[49][i] = 5; } for (k = 0; k < 10000; k++) for (j = 1; j < 49; j++) for (i = 1; i < 49; i++) { fc = 0.25 * (pot[i - 1][j] + pot[i + 1][j] + pot[i][j - 1] + pot[i][j + 1]); fs = 0.25 * (pot[i + 1][j + 1] + pot[i + 1][j - 1] + pot[i - 1][j - 1] + pot[i - 1][j + 1]); pot[i][j] = 0.8 * fc + 0.2 * fs; } for (i = 0; i < 50; i++) for (j = 0; j < 50; j++) printf("%d %d %f\n", i, j, pot[i][j]); return 0; } [/sourcecode]

Vertical Circular Motion with Friction

Introduction

This post analyzes the motion of a body in a vertical circular path, and the effect of friction therein. The calculations are performed for a value of g/r = 1, i.e., a radius of 10 m, taking the acceleration due to gravity to be 10 m/s². This should not materially alter results, except for the absolute value of time required to cover a certain angle.

Circular Motion without Gravity

Circular motion with friction without gravity is easily analyzed: the normal reaction is given by where r is the radius of the circular path, and v is the instantaneous velocity. We have which gives us where is the friction coefficient and is the angle through which it travels in the circle. Alas, there is no such happy conclusion for the motion where gravity is present. We are reduced to using numerical methods, some of which I outline below. The source code for a small test program is provided.

Frictionless Motion under Gravity

A body moving under gravity may be a pendulum that is released from a horizontal position, or a point mass released inside a smooth hemispherical bowl or cylinder. The conservation of energy gives the first integral, where we can get the position of the object in terms of the angle traversed θ. The elliptic integral formed can be solved to get the time required to reach any position. Note that a similar analysis can be made for a mass attached to the end of a light rod. However, the top point is a position of unstable equilibrium, so that the time taken to come down from that point in infinite.

Motion Under Gravity with Friction

Variation of Relative Time to Bottom with Friction

The above graph shows the variation of time required for a body released from the top of a circular path to reach the bottom. The times are normalized with respect to the time taken by a body to slide down a frictionless path.

Angular Velocity vs Angle for Vertical Circular Motion with and without Friction

The above graph shows the variation of angular velocity for a particle in vertical circular motion with and without friction. In the absence of friction, the angular velocity continuously increases, as the potential energy is converted into kinetic energy. Note that in the presence of friction, the maximum angular velocity need not occur at the bottom.

Decay of the Oscillation in Vertical Circular Motion with Friction

The above shows the motion of a body undergoing vertical circular motion with friction. The amplitude of the oscillation decreases each time. When the body is at the extreme position, if the angle made by the tangent to the circle at that path is less than the angle of repose, it will no longer slide down. Note that this means it is very difficult to make a body oscillate with small amplitude in the presence of friction. Also, the ideas on taking the friction force as approximately constant by considering small amplitudes don’t work either.

Where will the body stop? It will come to rest as the kinetic energy decreases, and this can occur at any point on the path. In fact, for the case under consideration, for a friction coefficient of more than 0.6, a body released from the horizontal position fails to reach the bottom.

It is possible to analytically solve the damped oscillation of a spring mass system, where the damping force is proportional to velocity. Whether the actual dampers found really obey this relation is a moot point. However, the characteristic equation can be analyzed to see whether the motion is overdamped, critically damped, or underdamped. The underdamped motion looks like the vertical circular oscillation plotted above.

Disclaimer: while the program is substantially correct, the code for determining the points where it reverses direction is not too good, so that there are some errors in the path.

#include
#include
#include

#define T 1.0

int main(int argc, char** argv){
float o, a, q;
float dt = T/100;
float mu;

int i;
int sgn;

mu = 0.1;
o = 0;
q = M_PI / 3;
sgn = -1;
for (i = 0; i < 1300; i++) { a = (cos(q) + sgn * mu * sin(q)) / T + sgn * mu * o * o; q += o * dt + 0.5 * a * dt * dt; o += a * dt; printf("%f %f\n", dt * i, -90 + q * 180 / M_PI); if (o < .02) { sgn *= -1; if (abs(cos(q)) < abs(mu * sin(q))) break; } } return 0; } [/sourcecode]

Projectile Motion with Drag

The projectile motion problem is the staple of high school physics classes, due mainly to its simplicity. However, the question arises, where do you go from here? How can you make sure that the calculated values and experiment agree? One approach is to model the other force acting on the body which is projected, the drag due to air.

Motion without Drag

Projectile motion for various angles with zero drag

The maximum range for a projectile for a given velocity of projection occurs at an angle of 45 degrees. The set of projectile paths for different angles at the same velocity will have an envelope that is a parabola.
The trajectory for a particular angle of projection θ is given by y = x tan θ – g x²/2u² cos ² θ, where u is the speed of projection.

Trajectory and Parameters of Motion

What is the relation between trajectory and the acceleration of the particle? We can have two particles with the same trajectory, but with entirely different accelerations. Consider this: is it enough for a particle to travel in a parabolic trajectory to simulate weightlessness? There are aircraft which do this, but they also have to match velocity. At the beginning of the climb, they should decrease speed while following the parabolic path till they reach the maximum height, and then increase the speed. A person following a parabolic path in a chopper moving slowly will not feel any difference in weight.

Motion with Drag

Projectile motion with high drag at different angles of projection

The drag acting on a projectile of a reasonable speed (say multiple of tens of metres per second) will be proportional to v^2, and the drag coefficient can be estimated based on some rough ideas on terminal velocity. With this data, the graphs for various angles are as shown. As can be seen, the range is much lower than that obtained from the drag free projection earlier.

Projectile motion with different drag coefficients

Consider the above case with a particle projected with different drag coefficients. The path, unlike the no drag case, is asymmetrical, with the motion till maximum height being at very high speeds, and the final path almost vertical. This is the typical path of a soccer ball. However, note that the actual drag on a body moving through air may not exactly be v². While it is comforting to have a simple characteristic like the drag coefficient describe the motion of a body, there are cases where it breaks down.

#include
#include
#include

int main(int argc, char** argv){
float u, q;
float ux, uy;
float dt = 0.05;
float vsq;
float x, y;
float ax, ay;
float coeff;
int i, j;

if (argc == 3) {
u = atof(argv[1]);
q = atof(argv[2]);
} else {
u = 100;
q = 45;
}

for (j = 0; j < 20; j++) { coeff = 0.016; q = 90 - j * 4; x = y = 0; printf("%f %f\n", x, y); ux = u * cos(q * M_PI /180); uy = u * sin(q * M_PI /180); ax = -coeff * u * u * cos(q * M_PI /180); ay = -coeff * u * u * sin(q * M_PI /180) - 10; for (i = 0; i < 1000; i++) { x += ux * dt; y += uy * dt; printf("%f %f\n", x, y); if (y < 0) break; ux += ax * dt; uy += ay * dt; vsq = ux * ux + uy * uy; ax = -coeff * vsq * cos(q * M_PI /180); ay = -coeff * vsq * sin(q * M_PI /180) - 10; } putchar('\n'); } return 0; } [/sourcecode]

Scores in Multiple Choice Papers

Score for 4:-1 scheme

Score 3:-1 scheme

#include
#include
#include
#include

#define NCH 4 /* Number of choices (4 or 5) */
#define LEN 30 /* Number of questions in the paper */
#define COR 3 /* Marks for correct */
#define WRG 1 /* Marks for wrong */

#define REPL 100000 /* Replications: set large enough */
#define BUCKETS 100

float urand(void);

int main(int argc, char **argv)
{
float rnd;
int i, j, k;

int soln[LEN]; /* The answer key */
int score; /* The marks scored */

int hist[BUCKETS];

for (i = 0; i < BUCKETS; i++) hist[i] = 0; /* Seed the random-number generator with current time so that * the numbers will be different every time we run. */ srand( (unsigned)time( NULL ) ); for (i = 0; i < LEN; i++) { rnd = NCH * urand(); /* U (0, NCH) */ for (j = NCH - 1; j >= 0; j–)
if (rnd > j) {
soln[i] = j;
break; /* Move to next rnd after bracket found */
}
}

for (k = 0; k < REPL; k++) { score = 0; for ( i = 0; i < LEN; i++) { rnd = NCH * urand(); /* U (0, NCH) */ for (j = NCH - 1; j >= 0; j–)
if (rnd > j) {
if (j == soln[i])
score += COR;
else
score -= WRG;
break;
}
}
hist[ BUCKETS * (score + LEN * WRG) /((COR + WRG) * LEN)]++;
}

for (i = 0; i < BUCKETS; i++) printf("%d %d\n", -LEN * WRG + (COR + WRG) * LEN * i / BUCKETS, hist[i]); putchar('\n'); return 0; } /* Uniform random number generator */ float urand(void) { return rand() * 1.0 / RAND_MAX; } [/sourcecode]

The above shows the distribution of marks obtained from a multiple choice answer scheme with 30 questions with 4:-1 and 3:-1 schemes. As can be seen, the 3:-1 is the correct scheme for such a paper, where a person guessing randomly will get zero credit. However, as can be seen, there is a small chance of getting positive scores. Certainly, it is not advisable to go into an exam with the idea of guessing all the answers. For an exam of reasonable difficulty, the chances of making it are practically zero. However, the more interesting question is this: if you have attempted most of the questions as well as you can, should you now guess the remaining question for a 3:-1 scheme. After all, on the average, you cannot lose, right? As the saying goes, in theory, theory and practice are the same.

Actual Test Data

Another interesting plot is as shown above: the actual data from a test looks very similar to the generated files in the previous section. What does this mean?

There is also an interesting application of the so called Monty Hall problem in the case of multiple choice questions. Consider a question with four options and only one correct answer. The probability of getting a correct answer from a random guess is 1/4. However, if you pick an answer, and later come to the conclusion that one of the other three is wrong, should you switch? The counter-intuitive answer is–yes! The probability of picking the correct answer if you switch is 1/2, rather than 1/3, if you had picked just one of the three options. That is to say, if you were given a three choice paper and picked one of them randomly, you would have 1/3 chance of getting it right. However, if you have four questions, in which you have eliminated one choice, then you have 1/2 chance of getting it right.

A	B	C	Keep	Switch
Dice	nil	nil	Win	Lose
nil	Dice	nil	Lose	Win
nil	nil	Dice	Lose	Win
			1/3	2/3

Let me elaborate–consider a test where you have four options A, B, C and D. Suppose you have picked option A. and then later found that one of the other options, say C is incorrect. Should you now switch?

The plots were generated using gnuplot.