ode45 goes to infinity!

Question

Lama Hamadeh 2020-3-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/509889-ode45-goes-to-infinity

评论： darova 2020-3-10

Dear all,

I am solving 2D shallow fluid equations where I first use the time-independent equations and then I integrate the time dependent to get the solution f or the vorticity. The problem is that ode45 goes to infinity and does not stop at the end of time domain. any ideas what might be wrong?

The main file is:

%Spatial variable on x direction
%--------------------------------
Lx=3; %spatial domain
delta=0.1; %spatial step size
xmin=-Lx; %minimum boundary
xmax=Lx; %maximum boundary 
Nx=(xmax-xmin)/delta; %number of spatial points
x=linspace(xmin,xmax,Nx); %spatial vector
%Spatial variable on y direction
%--------------------------------
Ly=3; %spatial domain
delta=0.1; %spatial step size
ymin=-Ly; %minimum boundary
ymax=Ly; %maximum boundary M
Ny=(ymax-ymin)/delta; %number of spatial points
y=linspace(ymin,ymax,Ny); %spatial vector
%Total matrix size
%-----------------
N = (Nx * Ny);
%Time variable
%-------------
dt=0.5;  %time step size
tmin=0; %minimum time
tmax=4; %maximum time
nt=(tmax-tmin)/dt; %number of time points
tspan=linspace(tmin,tmax,nt); %time variable
% base vector of zeros
e0 = zeros(N, 1); 
% base vector of ones
e1 = ones(N, 1);        
e2 = ones(N, 1); e2(mod(0:N-1, Nx) == 0) = 0; 
e3 = zeros(N, 1); e3(mod(1:N, Nx) == 0) = 1;
%%% Create matrix A %%%
ind_a = [1 Nx-1 Nx (N-Nx)];
diag_a = [e2 e3 e1 e1];
A = (1 / delta^2) * ...
    spdiags([rot90(diag_a, 2) -4*e1 diag_a], [-fliplr(ind_a) 0 ind_a], N, N);
%%% Create matrix B %%%
ind_b = [Nx N-Nx];
diag_b = [e1 -e1];
B = (1/(2*delta)) * ...
    spdiags([-1*rot90(diag_b, 2), diag_b], [-fliplr(ind_b), ind_b], N, N);
%%% Create matrix C %%%
e4 = repmat([0; ones(Nx-1, 1)], Nx, 1);
e5 = repmat([zeros(Nx-1, 1); -1], Nx, 1);
ind_c = [1, Nx-1];
diag_c = [e4, e5];
C = (1/(2*delta)) * ...
    spdiags([-1*rot90(diag_c, 2), diag_c], [-fliplr(ind_c), ind_c], N, N);
% Full matrices 
%--------------
A1 = full(A);
A2 = full(B);
A3 = full(C);
%Create a meshgrid to let Matlab know the x and y directions
%-----------------------------------------------------------
[X,Y] = meshgrid(x,y);
% Defining vorticity initial state/condition:
%--------------------------------------------
%2D Gaussian Function 
w0 = exp(-(X).^2 - ((Y.^2) / 20));
w_vec = reshape(w0, N, 1);
%solve for the vorticity
[Time,Omega]=ode45('vorticityBackslash',tspan,w_vec,[],A1,A2,A3); 

The function file is:

function rhs = vorticityBackslash(tspan,w_vec,dummy,A1,A2,A3) 
psi0 = A1\w_vec;
psix   = A2 * psi0;
psiy   = A3 * psi0;
wx     = A2 * w_vec  ;
wy     = A3 * w_vec  ;
diff_w = A1 * w_vec  ;
%Define viscosity
%----------------
visocity = 0.001;
%compute vorticity at a future step
%---------------------------------
rhs = visocity .* diff_w - psix.*wy + psiy.*wx;
end