Setting up a function to use with ODE solver

7 views (last 30 days)
Sean
Sean on 7 Mar 2013
Hi. I'm working on a project using one of the ode solvers in matlab to generate a numerical solution to an ode for charged particle motion. I've used the ODE solver before but I'm just confused on how to set up my function for my equation. The equation I'm using to solve for general motion for the charged particle is:
m*(dv/dt) = q(v x B)
I know that the general analytic solution to this problem is r = mV/qB but I am unsure of how to work this out with the ODE solution. If anyone could give me some help setting up my function it would be greatly appreciated. Thank you.
  6 Comments
Sean
Sean on 7 Mar 2013
Edited: Sean on 7 Mar 2013
This is what I have so far, I split up the velocity into x and y components:
function drdt = diffrv(t, rv)
q = 1.60217646 * 10^-19; % coulombs charge of electron
m = 9.10938188 * 10^-31; % mass of electron
B = 10^-9;
drdt = zeros(4,1);
drdt(1) = (q/m)*rv(2)*B;
drdt(2) = -(q/m)*rv(1)*B;
drdt(3) = (q/m)*drdt(2)*B;
drdt(4) = -(q/m)*drdt(1)*B;
drdt = [drdt(1);drdt(2);drdt(3);drdt(4)]
end
And for the actual solving:
init(1) = 5000; % x component velocity
init(2) = 3000; % y component velocity
init(3) = 0;
init(4) = 0;
tspan = 0:1:100;
options = odeset('RelTol', 1e-11, 'AbsTol', 1e-11);
[t,rdot] = ode45('diffrv',tspan,init,options)
plot(rdot(:,3),rdot(:,4))
end
Sean
Sean on 7 Mar 2013
This seems to work correctly also but it takes an insanely long amount of time to run

Sign in to comment.

Answers (2)

Babak
Babak on 7 Mar 2013
You need to write the equation in the dimentionless form. You cannot just simply write numbers like
m = 9.10938188 * 10^-31
in MATLAB and excpect it work fine!

Jan
Jan on 7 Mar 2013
Edited: Jan on 7 Mar 2013
How much faster is this:
function drdt = diffrv(t, rv)
q = 1.60217646e-19; % Avoid expensive power operation
m = 9.10938188e-31; % mass of electron
B = 1e-9;
c = (q/m)*B;
drdt = zeros(4,1);
drdt(1) = c*rv(2);
drdt(2) = -c*rv(1);
drdt(3) = c*drdt(2);
drdt(4) = -c*drdt(1);
% Useless: drdt = [drdt(1);drdt(2);drdt(3);drdt(4)]
end
When the relative and absolute tolerances are such tiny, long computing times can be expected. Unfortunately the results are not necessarily better, because small tolerances lead to a high number of integration steps and the rounding errors accumulate.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!