MATLAB: Discontinuous ODE . Why solvers dont reduce step size

discontinuousedo

Ive been tryng to solve an edo, but no matter which solver I use ( ode45, ode23, ode113, ode15s, ode23s,…) i ve got a wrong solution. The equation is:

dy = (Qin – Qout) / (pi * Rt ^ 2); Qout = siphon * C * sqrt(2 * g * y(1)) * pi * r ^ 2; Siphon can be 0 or 1. I ve got a solution using a solver with a fixed small step size. My question is: why the solvers were not able to reduce the step size when the variable siphon changes his value? thanks

Best Answer

See http://www.mathworks.com/matlabcentral/answers/59582#answer_72047. Matlab's integrators cannot handle non-smooth functions.

What a pitty that TMW included an example with a linear interpolation in the docs: https://www.mathworks.com/help/matlab/ref/ode45.html#bu3l43b . Not smart.

Related Solutions

MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

Here's a version:

% Initial conditions
t0=0;
x0=0:0.2:1;
y0=1:0.2:2;
dt=0.1;
tz=1;
t_range= t0:dt:tz;
X = zeros(numel(x0),numel(t_range));
Y = zeros(numel(y0),numel(t_range));
x_rk=zeros(1,numel(t_range));
y_rk= zeros(1,numel(t_range));
for n = 1:numel(x0)
    x_rk(1)=x0(n);
    y_rk(1)=y0(n);
    for i=1:length(t_range)-1
        [x_rk(i+1), y_rk(i+1)]=rk_method(x_rk(i),y_rk(i),dt);
    end
    X(n,:) = x_rk;
    Y(n,:) = y_rk;
    %% PLots
    subplot(numel(x0)/2,2,n)
    plot (t_range, x_rk,'b-o',t_range, y_rk,'y-o','MarkerSize',5);
    grid
    title([num2str(x0(n)),'   ',num2str(y0(n))])
    legend('x','y')
end
disp('X = '), disp(X)
disp('Y = '), disp(Y)
%% Functions
function [dxdt, dydt]= f(x,y)
dxdt =  5*x-3*y;
dydt = -6*x+2*y;
end 
function [x1, y1]=rk_method(x0,y0,dt)
[k11, k12]=f(x0,y0);     
[k21, k22]=f(x0+0.5*dt*k11,y0+0.5*dt*k12); 
[k31, k32]=f(x0+0.5*dt*k21,y0+0.5*dt*k22); 
[k41, k42]=f(x0+dt*k31, y0+dt*k32); 
x1=x0+dt*((k11/6)+((k21+k31)/3)+(k41/6));
y1=y0+dt*((k12/6)+((k22+k32)/3)+(k42/6));
end

MATLAB: Water Tank input that changes based on volume.

Please check with the following

clear, clc
V0 = 50;        %initial volume (m^3)
c0 = 0;     %initial concentration
Qin = 1.25 * 60;    %inflow (m^3/hour)
Qout = 1.1 * 60;    %outflow (m^3/hour)
cin = 100;          %pollutant concentration (mg/L)
t_final = 50;       %hours
dt = 0.5;           %time step
T = [dt:0.5:50];
t = 0;
%

k = 0;
%
%Initialize Loop Variables
cold = c0;
vold = V0;
inflow=1; % inflow status variable
for n = [1 : (t_final/dt)]
    t = t + dt;
    k = 0;
    
    if vold>=200 || inflow==0
        inflow=0; % say inflow stopped
        if vold<=150
            inflow=1;  % say inflow start
        end
    end
    
    if inflow==0
        Qin=0;
    else
        Qin=1.25 *60;
    end
    
    alpha = (Qin / (vold + (Qin - Qout) * t)) + k;
    beta = (Qin*cin)/(vold+(Qin - Qout) * t);
c(n) = cold - dt * alpha * cold + dt * beta;
V(n) = vold + (Qin - Qout) * dt;
 
cold = c(n);
vold = V(n);
end
yyaxis left
plot(T,V)
ylabel('Volume (m^3)')
xlabel('Time (hours)')
yyaxis right
plot(T,c)
ylabel('Concentration (mg/L)')
disp('      t        V           c')
format compact
disp([T' V' c'])
more(10)

Best Answer

Related Solutions

MATLAB: 4th order Runge-Kutta code that can solve for several intial conditions

MATLAB: Water Tank input that changes based on volume.

Related Question