Some weird results

[oysthaga]

I have a weird result for the relative error in problem 8, where the error have giant spikes at certain times: error.pdf error zoom.pdf I am pretty sure that I am evaluating the analytic and numeric values at the correct time.

And for the particles left, I get all the particles left for low frequencies, and none for high frequencies, even though the duration is much lower than 500 microseconds: Main.exe 7 1.25e-2 20, 100 particles.pdf (100 particles, timestep 1.25e-2 microseconds, duration 20 microseconds)

Main.exe 7 1.5625e-3 20, 10 particles.pdf (10 particles, timestep 1.5625e-3 microseconds, duration 20 microseconds)

7.8125e-4 20 , 10 particles.pdf (10 particles, timestep 7.8125e-4 microseconds, duration 20 microseconds)

The first of these took a couple of hours to run, so I can't use smaller timesteps and longer duration for 100 particles.

Any ideas what could be causing these issues?

[anderkve]

Hi @oysthaga!

Hmm, it's a bit difficult to guess just from the plots. If you share the link to your code on GitHub I might have time to take brief look to see if I spot anything suspicious.

One thing that might help you in finding the problem is to run a very short simulation, e.g. up to the first spike in the relative error plot, or even just a few time steps, and look at plots of x vs t, y vs t, and z vs t, with the graphs for the numerical and analytical solution in the same plot. That might help you spot if one of them seems to behave strangely.

Keep in mind that spikes in relative error can arise from a division by 0 (or close to 0). So what may be happening is that you have some actual bug that causes some generally larger-than-expected absolute error, and that this turns into huge spikes in relative error at the certain time steps when the analytical solution has coordinates close to 0.

[oysthaga]

Thought I linked the code, apparently I forgot. https://github.com/oysthaga/Project-3 The analytical values are defined in the Python file after line 76.

[oysthaga]

I just realized I forgot to include the time-dependent perturbation of the field in the analytical solution.

[anderkve]

You're not supposed to include it. The comparison with the analytical solution is in problem 8. We don't introduce the time-dependent E field until problem 9.

[oysthaga]

Right, I got confused for a minute, haha.

[anderkve]

I found a potential issue in computation of the relative error: In your Python code you have lines like this:

errEu0 = np.linalg.norm((r0 - r1Eu0)/r0, axis=0)

What this is doing is first dividing two 3-vectors r0-r1Eu0 and r0, which of course isn't mathematically defined, but numpy interprets this as element-wise division, so in the code you end up with a new 3-vector which contains the (signed) relative error per coordinate. And then you take the norm of that vector. However, the way we usually define the relative error for two vectors is as the ratio of two norms, e.g.

errEu0 = np.linalg.norm(r0 - r1Eu0, axis=0) / np.linalg.norm(r0, axis=0)

This can't explain the main issue, since we see that all the particles are escaping the trap with the time-dependent field, which is something unrelated to the computation of relative errors. But it may be worth looking at the relative error plots for the above definition of the relative error.

[oysthaga]

Oh, I see! Now I get this figure for the error:

New_error.pdf I think this seems reasonable?

[anderkve]

I think I found some issues with the RK4 code. You have code looking like this:

    std::vector<Particle> p = Particles; // Copies the particles because we need the old value in each step. 

    // [other stuff]

    for (int j = 0; j < p.size(); j++)
    {
        K1r[j] =  p[j].v*dt ;
        K1v[j] =  (total_force(j, t0, CoulombOn, omV, f)/p[j].m)*dt ;
        p[j].r = Particles[j].r + K1r[j]/2;
        p[j].v = Particles[j].v +  K1v[j]/2;
    }

There are two issues here, I think:

1) Since this is done in a single loop, the state of particle p[0] is updated before the K1 vectors are computed for particle p[1], and then the p[1] is updated before the K1 vectors are computed for particle p[2], etc. Since the computation of the K vectors depend on the state of the other particles, what you need to do is rather split this into two loops: First a loop that computes the K vectors for all the particles in their current state, and then a loop that updates the states of all particles using the K vectors. (The same of course applies for the loops of the subsequent steps of the RK4 algorithm.)

2) More importantly: In your loops you are updating the positions and velocities of the particle instances in the local vector p. But in computing K vetors you call the total_force function, which uses the current state of the particle instances living in the PenningTrap::Particle vector. That means that you are using the wrong positions and velocities. What you rather should do is, at each step in the RK4 algorithm, update the state of the particle in your PenningTrap object (i.e. the particles in Particles). So I think you should e.g. have lines like this:

        Particles[j].r = p[j].r + K1r[j]/2;
        Particles[j].v = p[j].v +  K1v[j]/2;

(Note that the first of these problems won't affect the relative error test, since in that case there is only a single particle.)

The new relative error plot indeed looks more sensible. Might I recommend a log-scale y axis? :)

Hope this helps -- I need to get back to other work now. Good luck finishing the project!

[oysthaga]

I will try this, thank you very much for the help. :)