Setting up a function to use with ODE solver

3 次查看(过去 30 天)
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 个评论
Sean
Sean 2013-3-7
编辑:Sean 2013-3-7
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 2013-3-7
This seems to work correctly also but it takes an insanely long amount of time to run

请先登录,再进行评论。

回答(2 个)

Babak
Babak 2013-3-7
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 2013-3-7
编辑:Jan 2013-3-7
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!

Translated by