ODE Solver Freezing Up

Question

Tom Keaton 2019-4-11

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/455840-ode-solver-freezing-up

编辑： Tom Keaton 2019-4-11

Hi,

I tried using ode15s to solve a trajecotry solution and what is so odd here is that the solver exponentially becomes slower as the time steps persist. I tried using the additional options to set the AbsTol and RelTol to a certain value and that did not help either. Can someone help me understand why the solver is freezing up in the middle? I used the tic-toc function to see whats going on and this was the output for the ode solver:

As you can see it gets exponentially slower when solving. This is the code:

%TSTEP MUST BE ABLE TO GO INTO TFIN AN ODD AMOUNT OF TIMES (Because simulation includes tstep @ 0 so it ends up being even)
tstep = 0.25E-8; %Defining time step
tfin = 10.25E-7; %Defining final time
ienergy = 1E-20;
%Initial conditions to solve for pitch angle and lambda value
lambdacheck = 30; %In degrees
%options = odeset('AbsTol',1E-7);
%Initial velocity and position
sine2alpha = (cos(lambdacheck)^6)/sqrt(1+3*sin(lambdacheck)^2);
m_e = 9.11E-31;
v0mag = sqrt(2*ienergy/m_e);
vx = 0;
vy = sqrt(sine2alpha)*v0mag;
vz = sqrt(1 - sine2alpha)*v0mag;
x = 0.05;
y = 0;
z = 0;
tic %Begin timer
%Natural Constants
m_e = 9.11E-31; %Mass of electron
jtoev = 6.24E+18; %Joules to eV conversion
epnaut = 8.854187E-12;
k = 1/(4*pi*epnaut);
%Create zeros matrices to populate later
energymat = zeros(particlecount,1);
vmagmat = zeros(particlecount,1);
intspan = [0:tstep:tfin]; %Total time span
intspancol = intspan(:); %inverts matrix's rows and columns
[introw,intcol] = size(intspancol);
icposition = zeros(1,3);
%Generate matrix that will be populated by positions and velocities (Each matrix within structure "Wposandvel" should be of size 2*introw-1 by 6)
Wposandvel = zeros(2*introw-1,6);
%This section is purely for generating matrices the program will populate later
tindex = [0:1:introw-2];
tindexcol = tindex(:); %inverts matrix's rows and columns
[tinrow,tincol] = size(tindex);
PCPmat = zeros(tincol,1); %Probability of collision solutions
choicemat = zeros(tincol,1); %Set of RNG choices system makes
colnewv = zeros(intcol,3); %New velocity components after collision
%Vsweepmat = zeros(intcol,1); %Volume swept out by particle in given time step
vmagmat = zeros(intcol,1); %Magnitude of velocity at given time step
evmat = zeros(intcol,1); %Energy of eV at every time step
%crossareamat = zeros(intcol,1); %Cross-sectional area at every time step
for t = 0:1:introw-2 %complete time interval and time step
        
      [trow,tcol] = size(t);
      r = sqrt(x.^2 + y.^2 + z.^2);
	vmag = sqrt(vx.^2 + vy.^2 + vz.^2);
	%Coupled differential equation solver for trajectory within current time step
	icv = [x; y; z; vx; vy; vz]; %Initial conditions
      tspan = [intspan(t+1) ((intspan(t+2)-intspan(t+1))/2)+intspan(t+1) intspan(t+2)]; %Time span
        
      tic
      %Trajectory solver
    	[T,S] = ode15s(@bdipuniodefun, tspan, icv);
    	[rownum,colnum] = size(S);
       Wposandvel((1+2*t):(3+2*t),(1:6)) = S;
       toc
        
	%Redfine the velocity and position components to reference on next if-loop run
      x = Wposandvel((3+2*t),1);
	y = Wposandvel((3+2*t),2);
	z = Wposandvel((3+2*t),3);
	vx = Wposandvel((3+2*t),4);
	vy = Wposandvel((3+2*t),5);
	vz = Wposandvel((3+2*t),6);
        
end

Here is the ode set the solver uses:

function bdip = bdipuniodefun(t,s)
%Using SI units
q = -1.60217662E-19;
m_e = 9.11E-31;
persistent Bx By Bz Ex Ey Ez
if isempty(Bx)
    [Bx, By, Bz] = B_test();
end
bdip = [s(4); s(5); s(6); (q/m_e)*(Ex(s(1),s(2),s(3)) + s(5)*Bz(s(1),s(2),s(3)) - s(6)*By(s(1),s(2),s(3))); (q/m_e)*(Ey(s(1),s(2),s(3)) + s(6)*Bx(s(1),s(2),s(3)) - s(4)*Bz(s(1),s(2),s(3))); (q/m_e)*(Ez(s(1),s(2),s(3)) + s(4)*By(s(1),s(2),s(3)) - s(5)*Bx(s(1),s(2),s(3)))];
end

Here is B_test:

function [Bx, By, Bz] = B_test()
Bfieldstrength = 0.64; %In (Teslas)
magvol = 3.218E-6; %In (m)
mu0 = (4*pi)*10^-7;
magnetization = (Bfieldstrength*magvol)/mu0;
syms x y z
m = [0,0,magnetization];
r = [x, y, z];
B = mu0 *(((dot(m,r)*r*3)/norm(r)^5) - m/norm(r)^3);
Bx = matlabFunction(B(1));
By = matlabFunction(B(2));
Bz = matlabFunction(B(3));
end