Setting up a function to use with ODE solver

6 ビュー (過去 30 日間)
Sean
Sean 2013 年 3 月 7 日
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.

カテゴリ

Help Center および File ExchangeNumeric Solvers についてさらに検索

Community Treasure Hunt

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

Start Hunting!

Translated by