ode45 goes to infinity!

6 views (last 30 days)
Lama Hamadeh
Lama Hamadeh on 9 Mar 2020
Commented: darova on 10 Mar 2020
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
  1 Comment
darova
darova on 10 Mar 2020
I don't why it's working. You are passing parameters in a wrong way
Here is how it should be
[Time,Omega]=ode45(@(t,w)vorticityBackslash(t,w,[],A1,A2,A3),tspan,w_vec);

Sign in to comment.

Answers (0)

Community Treasure Hunt

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

Start Hunting!