*Concern highlighted in yellow

*Solution at bottom

I have a differential equation to solve for the motion of an electron:

$$ \frac{d^2v}{dt^2} = \frac{1}{\gamma^6}\left( \frac{eE}{\tau m} - \left( \frac{\gamma}{\tau} \frac{dv}{dt} + \frac{\gamma^3}{\tau c^2} v^2 \frac{dv}{dt} + 6\frac{\gamma^8}{c^2} v \left( \frac{dv}{dt} \right)^2\right) \right) $$

$$ \gamma = \frac{1}{\sqrt{1-v^2/c^2}}\\ e = 1.6\times 10^{-19}\\ E = 10^6\\ \tau = 6\times 10^{-24}\\ m = 9.1\times 10^{-31}\\ c = 3\times10^8 $$

With initial conditions $$v(0) = 0$$ $$\frac{dv}{dt}\biggr\rvert_0 = Ee/m \approx 6\times 10^{-4}$$

But I heard that it is good practice to make everything dimensionless so that your computer can handle things nicely, so I introduced some dimensionless variables:

$$ T = t/\tau\\ V = v/c\\ \varepsilon = \frac{\tau e}{cm} E $$

This gets me to the final form:

$$ \frac{d^2V}{dT^2} = \frac{1}{\gamma^6}\left( \varepsilon - \left( \gamma \frac{dV}{dT} + \gamma^3 V^2 \frac{dV}{dT} + 6\gamma^8 V \left( \frac{dV}{dT} \right) ^2 \right)\right) $$

With initial conditions

$$ V(0) = 0\\ \frac{dV}{dT}\biggr\rvert_0 = \varepsilon \approx 4\times 10^{-15} $$

I am using MATLAB's ode45 to solve this, but I do not think it can handle such small initial conditions. The plot should rapidly asymptote to 1, but I can only get that to happen with fairly large initial conditions:

Initial acceleration equal to $\varepsilon$

enter image description here

Initial acceleration equal to $10$

enter image description here

I can't tell what the issue is, but I'll post my code - it's not very long:

Script for the Differential Equation

function dv2dt = emDeq(t, v)
mu = 4*pi*10^(-7);
c = physconst('LightSpeed');
e = 1.6022e-19;
m = 9.1094e-31;
tau = mu*e^2/(6*pi*m*c);

E = (10^6)*(tau*e)/(c*m);
% E = 1e6;

gam = (1 - v(1)^2)^(-1/2);

dv2dt = [v(2); (gam^(-6))*( E  - ( gam*v(2) + (gam^3)*(v(1))^2*v(2) + 6*(gam^8)*v(1)*(v(2))^2 ) )];


Script for Plotting

mu = 4*pi*10^(-7);
c = physconst('LightSpeed');
e = 1.6022e-19;
m = 9.1094e-31;
tau = (mu*e^2)/(6*pi*m*c);
E = (10^6)*(tau*e)/(c*m);
% E = 1e6;

v0 = 0;
a0 = E;
% a0 = E*e/m;
[t, v] = ode45(@emDeq, [0, 10], [v0, a0]);
figure('Color', 'w')
plot(t, v(:,1));
title('V vs T'); xlabel('T'); ylabel('V');


I needed to redefine length scales and time scales, as well as use either ode15s or ode23s.

1 time step is $\tau$s, one length step is $10^{-15}$m, one voltage step is $10^{-10}$. This got me a reasonable $\epsilon$ to work with. Plots below. Thanks to everyone who helped out!

enter image description here

  • $\begingroup$ You desperately need to either change units or nondimensionalize your problem. The magnitudes on the constants and dependent variable are likely causing trouble. I would also recommend you use a stiff solver like ode15s just to be sure that these constants are not causing stability issues that can be easily solved $\endgroup$ – whpowell96 Dec 2 '19 at 19:19
  • $\begingroup$ @whpowell96 I nondimensionalized it, but the magnitude issues seemed to have gotten compressed to a single constant, $\varepsilon$, which appears in the initial conditions as well as the differential equation. I am not sure how to deal with this. $\endgroup$ – MurderOfCrows Dec 3 '19 at 3:17
  • $\begingroup$ Maybe you should try rescaling the time like $s=\varepsilon^\alpha T$, $U = \varepsilon^\beta V$. Rewrite the equation in terms of $U$ and $s$ then you can play around with $\alpha$ and $\beta$ to fix to rescale the initial condition to your liking. It may also be possible to choose $\alpha$ and $\beta$ so it makes the actual differential equation better conditioned as well. Again, I would suggest using ode15s just in case $\endgroup$ – whpowell96 Dec 3 '19 at 4:16
  • $\begingroup$ I worked with your suggestion and I got somewhere hopeful! I will post the resolved version. $\endgroup$ – MurderOfCrows Dec 3 '19 at 4:37
  • $\begingroup$ Great! Don't forget to rescale everything back in the end to make sure your units are correct $\endgroup$ – whpowell96 Dec 3 '19 at 5:01

I think the problem with initial conditions is more related with the sum of very small numbers with other that are only small (I suspect the problem resides in terms like gam^8, since gam is small). The resulting numeric inaccuracy will bring problems.

I've played around with other initial conditions and integration times and in octave (that is similar but not equal to matlab), the integration seems to work initially for $E=10^8$ but not for $E=10^7$, although I observed numeric instabilities for $t>10$. That numeric instability was observed up to $E=10^{16}$ but not for $E=10^{17}$.

Although it does not solve your problem in this case, sometimes changing the order of the parcels may help in reducing the numeric inaccuracy problems, e.g., something like:

dv2dt = [v(2); -(gam^(-6))*( 6*(gam^8)*v(1)*(v(2))^2 + (gam^3)*(v(1))^2*v(2) + gam*v(2) - E )];
  • $\begingroup$ gam is at least 1 and blows up fairly quickly, so I am not sure if that is the issue. Though relative sizes may be, so perhaps I need to use a more precise ode solver. Do you know of any? $\endgroup$ – MurderOfCrows Dec 2 '19 at 15:44
  • $\begingroup$ @MurderOfCrows I'm not sure if your problem is stiff or not, but perhaps ode15s. $\endgroup$ – Ertxiem - reinstate Monica Dec 2 '19 at 23:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.