3

こんにちは、さまざまな力でy方向に移動するx平面のforループを使用して、さまざまな粒子を作成できる小さなプログラムを書いています。粒子は反発力で相互に作用するため、ode45 の計算では、粒子相互の位置が必要になります。コードを何度も見直しましたが、解決策が "NaN" として表示されるため、どこに問題があるのか​​ わかりません。時間間隔は、すべての粒子の位置が ode45 計算で更新されるように、粒子の動きが時間内に分割されるようになっています。時間を分割するためにイベントを使用しました。コードがあまり明確でないため、さらに情報が必要な場合は、お問い合わせください。コードに定数を示していません。

time_interval = 1000; %Time in microseconds
%The number of time intervals the particle movement is divided by
number_of_intervals = 10000/time_interval;

tstart=0;
tstop= time_interval*1e-6;

num_of_particles = 3;

w0{1} = [0;0;0;0];
sols{1} = transpose(w0{1});
x_adj_stable(1) = 0;
y_adj_stable(1) = 0;

break_vect(num_of_particles) = 0;
times{1} = 0;

%Initiate particle characteristics
for particle_num = 2:num_of_particles
w0{particle_num} = [0;0;w0{particle_num - 1}(3) + 4*1e-6;0];
x_adj_stable(particle_num) = w0{particle_num}(3);
y_adj_stable(particle_num) = 0;
sols{particle_num} = transpose(w0{particle_num});
times{particle_num} = 0;
end

%Event options, Identifying the time intervals within the solution
options = odeset('Events', @events);

%The loop which connects the intervals of time
for n = 1 : number_of_intervals

    for particle_num = 1:num_of_particles
%===========================================
[t_sol, y_sol] = ode45(@pm1dwoc, [tstart, tstop], w0{particle_num}, options);
%===========================================
%The number of steps within each interval
n_steps = length(t_sol); 

for interval = 2: n_steps
        times{particle_num} = [times{particle_num}; t_sol(interval)];
        sols{particle_num} = [sols{particle_num}; y_sol(interval, :)];
end

%Setting the initial conditions of the next interval

w0{particle_num} = y_sol(n_steps,:);

y_adj(particle_num) = y_sol(n_steps, 1);

x_adj(particle_num) = y_sol(n_steps, 3);

    end

tstart = tstart + time_interval*1e-6;
tstop = tstart + time_interval*1e-6;

for n = 1:num_of_particles
y_adj_stable(n) = y_adj(n);
x_adj_stable(n) = x_adj(n);
end
end

 function dwdt = pm1dwoc (t,w)
    y=w(1);
    vy=w(2);
    x=w(3);
    vx=w(4);

    for particle_number = 1:num_of_particles 
    Dist(particle_number) = sqrt((x - x_adj_stable(particle_number))^2 + (y - y_adj_stable(particle_number))^2);

    if (Dist(particle_number) > 0)
        x_vect(particle_number) = (x - x_adj_stable(particle_number))/Dist(particle_number);
        y_vect(particle_number) = (y - y_adj_stable(particle_number))/Dist(particle_number);
    end
    end


     dy_component = -qw*Vd/(m*G) + qw*(-qw)/(m*4*pi*epsilon0*(2*R+2*y)^2));
     dx_component = 0;
     for particle_number = 1:num_of_particles
         if (Dist(particle_number) > 0)
            dx_component = dx_component + ((qw^2)/(0.0001*m*4*pi*epsilon0*Dist(particle_number)^2))*x_vect(particle_number);
            dy_component = dy_component + ((qw^2)/(0.0001*m*4*pi*epsilon0*Dist(particle_number)^2))*y_vect(particle_number);
         end
     end

     dwdt = [vy; dy_component; vx; dx_component];

end

function [eventvalue,stopthecalc,eventdirection] = events (t,w)
    t= t*1e6;
    end_time = 10000;

    eventvalue = t - end_time/number_of_intervals;
    stopthecalc = 1;
    eventdirection = 1;

    for num = 2:number_of_intervals
        eventvalue = [eventvalue; t - num*(end_time/number_of_intervals)];
        stopthecalc = [stopthecalc; 1];
        eventdirection = [eventdirection; 1];
    end

end

end

どんな助けでも本当に感謝します!

ありがとうございました

サム

編集

バグがないことがわかっているセグメントを取り出したので、簡略化したバージョンを投稿します。

num_of_particles = 3;

w0{1} = [0;0;0;0];
sols{1} = transpose(w0{1});
x_adj_stable(1) = 0;
y_adj_stable(1) = 0;
times{1} = 0;

%Initiate particle characteristics
for particle_num = 2:num_of_particles
    w0{particle_num} = [0;0;w0{particle_num - 1}(3) + 4*1e-6;0];
    x_adj_stable(particle_num) = w0{particle_num}(3);
    y_adj_stable(particle_num) = 0;
    sols{particle_num} = transpose(w0{particle_num});
    times{particle_num} = 0;
end

for particle_num = 1:num_of_particles

    [t_sol, y_sol] = ode45(@pm1dwoc, [tstart, tstop], w0{particle_num});

    %The number of steps within each interval
    n_steps = length(t_sol); 

    for interval = 2: n_steps
            times{particle_num} = [times{particle_num}; t_sol(interval)];
            sols{particle_num} = [sols{particle_num}; y_sol(interval, :)];
    end

    %Setting the initial conditions of the next interval
    w0{particle_num} = y_sol(n_steps,:);
    y_adj(particle_num) = y_sol(n_steps, 1);
    x_adj(particle_num) = y_sol(n_steps, 3);

end

%//This part is written since the y_adj and x_adj variables are constantly changing so              %//stable variable version is required when used in the ode45 calculations
for n = 1:num_of_particles
    y_adj_stable(n) = y_adj(n);
    x_adj_stable(n) = x_adj(n);
end
end

 function dwdt = pm1dwoc (t,w)
    y=w(1);
    vy=w(2);
    x=w(3);
    vx=w(4);

    for particle_number = 1:num_of_particles 
    Dist(particle_number) = sqrt((x - x_adj_stable(particle_number))^2 + (y - y_adj_stable(particle_number))^2);

    if (Dist(particle_number) > 0)
        x_vect(particle_number) = (x - x_adj_stable(particle_number))/Dist(particle_number);
        y_vect(particle_number) = (y - y_adj_stable(particle_number))/Dist(particle_number);
    end
    end


     dy_component = 0;
     dx_component = 0;
     for particle_number = 1:num_of_particles
         if (Dist(particle_number) > 0)
            dx_component = dx_component + (1/(Dist(particle_number))*x_vect(particle_number);
            dy_component = dy_component + (1/(Dist(particle_number))*y_vect(particle_number);
         end
     end

     dwdt = [vy; dy_component; vx; dx_component];

end
4

1 に答える 1

1

一般に、特に方程式が連成している場合は、多くの方程式を個別に積分することは望ましくありません。ode45そしてその親戚は、連立方程式系を完全に統合することができます。

簡単な例として、3 つの変数の導関数を定義する次の関数があるとします。派生物は互いに依存しています。

function xprime = my_ode(t, x)

xprime = [x(1) + 2*x(3);
          x(2) - x(1);
          3*x(3) + x(1) + x(2)];

これはあまり興味深いシステムを生成しませんが、3 つの異なる変数を一度に統合する方法を示しています。を使用するode45と、次のように入力できます。

x0 = [0;0;1];
tspan = 0:0.1:1;
[T,X] = ode45(@my_ode, tspan, x0)

問題の場合、この構造により、各ステップで粒子間の分離を計算し、反応を決定できます。次に、それぞれを統合して効果を調整する代わりに、それらの効果を運動方程式の一部にすることができます。

于 2012-12-05T07:17:16.410 に答える