Solving vectorized ODE (Solving same ODE with many initial conditions at once).
显示 更早的评论
I am tring to apply 2nd newton law to many points, where m,r,v,f is a n-by-1, n-by-3, n-by-3, n-by-3, n-by-3, matrixes for mass, position,velocity, and force for n points. The 3 columns are the x,y,z components for those values.
options = odeset('Vectorized','on');
y0=[r,v];
[t,YSol]=ode45(@(t,y) MotionODE(t,y,m,f),[0,dt],y0,options);
r=YSol(t==dt,1:3)';
v=YSol(t==dt,4:5)';
function d2rdt2= MotionODE(t,Y,m,f)
%ODE function for motion
r=Y(:,1:3);
v=Y(:,4:5);
drdt=v;
dvdt=f./m;
d2rdt2=[drdt,dvdt];
end
However I got an error: "Index in position 2 exceeds array bounds (must not exceed 1).". Is the option "vectorized" designed for what I think it is for? How do I get this run correctly (other than a for loop)?
采纳的回答
Walter Roberson
2021-3-17
r = Y(1,:);
v = Y(2,:);
drdt = v;
dvdt = f./m * ones(size(v));
d2rdt2=[drdt;dvdt];
3 个评论
Yi-xiao Liu
2021-3-17
Walter Roberson
2021-3-17
Okay, I misread earlier. 'Vectorized' is not useful for your situation.
%m is n x 1
%r is n x 3
%v is n x 3
%f is n x 3
%the boundary conditions are arranged in memory as
%all of the position x coordinates, then all of the position y, then all
%of the position z, then all of the velocity x, then all of the velocity y,
%then all of the velocity z
y0 = [r,v]; %[n x 3, n x 3] --> n x 6
f_over_m = f ./ m; %n x 3 ./ n x 1 -> n x 3
[t,YSol] = ode45(@(t,y) MotionODE(t, y, f_over_m), [0,dt], y0);
R = reshape(YSol(end,1:end/2), [], 3); %px, py, pz
V = reshape(YSol(end,end/2+1:end), [], 3); %vx, vy, vz
function d2rdt2= MotionODE(t, Y, f_over_m)
%ODE function for motion
Y = reshape(Y, [], 6); %px, py, pz, vx, vy, vz
r = Y(:,1:3); %n x 3
v = Y(:,4:5); %n x 3
drdt = v; %n x 3
dvdt = f_over_m; %n x 3
d2rdt2 = reshape([drdt, dvdt], [], 1); %MUST return column vector
end
Yi-xiao Liu
2021-3-17
更多回答（1 个）
Jan
2021-3-17
1 个投票
same ODE with many initial conditions at once
This is a bad idea. Remember that the step size control is triggered by the most susceptible component of the trajectory. This reduces the stepsize for all components and in consequence increases the accumulated rounding errors. The total number of calculations can be larger than running the integration for each initial value separately.
4 个评论
Yi-xiao Liu
2021-3-17
Jan
2021-3-18
The error can be small, if the trajectory is not very sensitive for a variation of the initial values. If you integrate a chaotic function like a double pendulum, the error will get much larger.
You can improve the performance by choosing another integrator or my modifying the tolerances. In your case the function to be uintegrated is tiny. Therefore saving function calls is not the standard way to improve the speed. If this function contains thousands of lines of code, the situation would be different.
Yi-xiao Liu
2021-3-18
Jan
2021-3-18
Do you have the Parallel Processing Toolbox? Then try to replace the FOR by PARFOR.
类别
在 帮助中心 和 File Exchange 中查找有关 Programming 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
