Please help solve a system of differential equation

I have a system of differential equation in x,y with boundary conditions. I want to solve it analytically and numerically. The analytical solutions coming in the form of error functions. Can someone please help me find analytical and numerical solution for this, I have been trying for weeks now, really need to solve this for my thesis. Please help.

Answer 1

[Alan Stevens]

Here's a numerical approach:

First manipulate the equations

tspan = [0, 120];
ic = [80000, 0.001];
[t, xy] = ode45(@rate, tspan, ic);
x = xy(:,1); y = xy(:,2);

subplot(2,1,1)
plot(t,x),grid
xlabel('t'), ylabel('x')
subplot(2,1,2)
plot(t,y),grid
xlabel('t'), ylabel('y')

function dxydt = rate(~,xy)

        a = 0.1; b = 0.2; c = -0.1;  % Replace with desired values
        
        x = xy(1); y = xy(2);
        dxydt = [a-b+c*x*y;            % eqn (3) 
                 c*y*(1-y) - a*y/x];   % eqn (5)
end

[Jayesh Jawandhia]

Hey Alan, thank you so much for the help. I was just hoping if there is something we can do to find the right solution, x and y will always be positive and increase with time initially. I also have the experimental data for x,y to validate. x and y both should ideally increase and reach a saturation. If you want i can also share the experimental data with you. It'll be really helpful if you could help with this, I have been trying to solve this for weeks.

[Alan Stevens]

Do you know the values of a, b and c, or do you need to fit them using the experimental data? If the latter you will need to upload the data if you want any help.

[Jayesh Jawandhia]

Hey Alan, a and b need to be fitted using the experimental data. c is not known can be taken according to the our requirement. Attaching the data for your reference. Thank you so much for all the help.

[Alan Stevens]

I'm puzzled! Looking at your data I see the initial value of x is likely to be around 10 to 15, not 80000. Where does the latter come from? The data file also lists different values of a and b for every data point. Are these just for guidance to typical values?

[Alan Stevens]

This data looks more like it's related to the logistic function than the error function. The following shows a purely empirical fit:

t = 5:5:120;
x = [15.6700   18.0000   19.0000   19.3300   19.5000   19.6700   19.3300 ...
     18.7700   18.6000   18.2000   17.9300   17.9300   17.6700   17.3300 ...
     17.1700   16.8300   16.8300   16.1700   15.5000   15.5000   15.8300 ...
     15.6700   15.6700   15.6700];
 
 y = [53.4467   61.4233   63.9467   64.6933   65.6633   65.7567   65.6100 ...
     65.3000   65.6133   65.9367   66.2833   66.0300    67.0500   65.9600 ...
     66.3733   66.7867   66.5233   66.3533   66.6567    66.6400   66.4267 ...
     67.3067   67.0300   67.3133];
 
 a = 67; b = 0.2; c = 0.01; abc = [a, b, c];
 abc = fminsearch(@(abc) fny(abc, t, y), abc);
 a = abc(1); b = abc(2); c = abc(3);
 Y = a.*exp(-c*t)./(1+exp(-b*t));
 
 A = 20; B = 0.2; C = 0.01; ABC = [A, B, C];
 ABC = fminsearch(@(ABC) fnx(ABC, t, x), ABC);
 A = ABC(1);B = ABC(2); C = ABC(3);
 X = A.*exp(C*t)./(1+exp(-B*t));
 
 plot(t,x,'ko',t,y,'rs',t,Y,'k',t,X), grid
 xlabel('t'),ylabel('x and y')
 legend('x','y','yfit','xfit')
 axis([0 120 0 90])
 hold on

 
 
 function Fy = fny(abc, t, y)
 
            a = abc(1); b = abc(2); c = abc(3);
            Y = a.*exp(-c*t)./(1+exp(-b*t));
            
            d = Y-y;
            Fy = sum(norm(d));
 end

 function Fx = fnx(ABC, t, x)
           
          A = ABC(1); B = ABC(2); C = ABC(3);
          X = A.*exp(C*t)./(1+exp(-B*t));
          d = X - x;
          Fx = sum(norm(d));
 
 end

[Jayesh Jawandhia]

Hey Alan, really sorry. I missed the units in the data. The initial value in the question was in a different unit. It should 12.5. I am attaching the correct data here. Also, a and b ideally will be constant so if we want we can take an average of them and use it in the numerical solution. During experiments, it was coming time dependent due to experimental errors. We want to get a numerical solution through equations. Thank you for all the help, it means alot. Please check this whenever you get some time.

[Alan Stevens]

1. What is the system you are trying to calculate? i.e. What do x and y physically represent?
2. How did you derive the equations, and why does your data set contain estimates for a and b, but not c?
3. How did you derive the point estimates for a and b in the data set?

I ask because your equations look strange and don't seem to represent the data, whereas your x and y data just look something close to logistic equations.

[Jayesh Jawandhia]

Hey Alan,

Trying to answer your questions below:

1. The system is trying to solve for the Volume of liquid produced in a reactor (x) and % of air present in that liquid (y). Basically the system is a continuous-stirred tank reactor with inlet and outlet feed. The inlet is reacting with air and a solution is formed with some air being trapped with the inlet.
2. We derived these equations using mass balance (Inlet - Outlet + Generation = Accumulation) on both air and liquid.
3. a and b are the inlet and outlet feed rates. These are being pumpd using a peristaltic pump but due to problems these are coming as a function of time. c is the rate constant in the interaction of air with the inlet to form the liquid.

Thank you so much for all the time, really appreciate it. Please let me know if this helps, happy to answer any other questions. I haven't been able to figure out why this is going wrong.

[Alan Stevens]

Hmm. Well, the approach I'd take is as follows:

1. Estimate initial values for a, b and c.
2. Numerically solve the equations to get predicted values at the experimental timesteps.
3. Form residuals at those timesteps: Experimental - predicted.
4. Calculate the sum of squared residuals.
5. Use fminsearch to try to find values of a, b and c that minimises the sum of squared residuals.

[Jayesh Jawandhia]

Hey Alan,

Thank you for the approach. I was hoping if you could help with the scenario when we take constant a and b. We are currently trying to solve for a constant a and b and get to the required curve behaviours.

[Alan Stevens]

I did have a quick go at this. Unfortunately, I wasn't able to make a sufficiently good estimate of intial values to get converged results!

[Jayesh Jawandhia]

Hey Alan, can you share that if possible? I'll try to work on that.

[Sam Chak]

@Jayesh Jawandhia, If the process system is controllable, can you inject a Manipulated Variable u so that you can freely change to see how the State Variables respond?

It will be way more flexible and practical to design (that drives the system to behave as desired) than to repeatedly guess a unique set of initial values to arrive at the desired converged result.

A Continuous Stirred Tank Reactor (CSTR) is a controllable process system.

[Jayesh Jawandhia]

Hey @Sam Chak,

Thank you for your response. I think we can do that, we can try to do that. I am not so familiar with these numerical methods, can you share some reference or material that I can try to use and solve it? It'll be really helpful. Really appreaciate it.

[Sam Chak]

I'm no expert in CSTR. I read your comments and Google the keywords. Since CSTR is a vessel in the industrial process, it is definitely controllable. Compare the equations in the journal papers and identify the manipulated variables, which can be more than 1. The math can come later.

More importantly, get familiar with the CSTR differential equations with the manipulated variables.

[Alan Stevens]

@Jayesh Jawandhia Ok. Here is the code I put together to have a quick look. Change the values of a, b and c to your heart's content! I'm not convinced that your odes describe your data, but more than happy to be proved wrong. If you are successful, post your results here. Good luck..

t = 0:300:7200;

x = [15.0000   17.2300   18.1900   18.5000   18.6600   18.8200   18.5000 ...
     17.9600   17.8000   17.4200   17.1600   17.1600   16.9100   16.5900 ...
     16.4300   16.1100   16.1100   15.4700   14.8400   14.8400   15.1500 ...
     15.0000   15.0000   15.0000];
 
 y = [0.5345    0.6142    0.6395    0.6469    0.6566    0.6576    0.6561 ...
      0.6530    0.6561    0.6594    0.6628    0.6603    0.6705    0.6596 ...
      0.6637    0.6679    0.6652    0.6635    0.6666    0.6664    0.6643 ...
      0.6731    0.6703    0.6731];
 
 % Add initial values to data
 x = [12.5 x];
 y = [0.001 y];
 
 % Initial guesses for a, b and c
 %      [a, b, c]
  abc = [1, 1, 1];
 
  % Try fminsearch to estimate better values of a, b and c
  ABC = fminsearch(@(abc) compare(abc, t, x, y), abc);
  
  disp(ABC)  % Display resulting values
  
  % Numerically integrate the odes with resulting values
  [t, XY] = runode(ABC,t);
  
  % Extract fitted values and plot fits vs data
  X = XY(:,1);  Y = XY(:,2);
  subplot(2,1,1)
  plot(t, X, 'r', t, x, 'ro'),grid
  ylabel('x')
  axis([0 7200 0 20])
  subplot(2,1,2)
  plot(t, Y, 'r', t, y, 'ro'),grid
  ylabel('y')
  axis([0 7200 0 1])
 
  % Calculate sum of squared residuals
  function F = compare(abc, t, x, y)
  
           [~, XY] = runode(abc, t);
           X = XY(:,1); Y = XY(:,2);
           dx = x - X;
           dy = y - Y;
           F = norm(dx)+norm(dy);
           
  end
  
  % Call the odes with current values of a, b and c
  function [t, XY] = runode(abc, t)
  
  ic = [12.5, 0.001];
  [t, XY] = ode45(@(t,xy) odefn(t, xy, abc), t, ic);
  end
  
% ODE equations
 function dxydt = odefn(~,xy, abc)
           a = abc(1); b = abc(2); c = abc(3);
           x = xy(1); y = xy(2);
           dxydt = [a-b+c*x*y;
                    c*y*(1-y)-a*y/x];
 
 end

[Jayesh Jawandhia]

Thanks @Alan Stevens @Sam Chak I'll try to put these into work and update you guys. Thank you so much, this is really helpful.

[Jayesh Jawandhia]

Hey @Alan Stevens

Hope you are doing well. I have been reading up on these and your solution has been really helpful. I just had a follow-up. We are trying to find a better fit by making a and b as a function of t. I haven't been able to figure out how do I add these to my ode45 solver in Matlab just like what you did here.

Can you please help with this? For simplicity, I am assuming both a and b second order in t, i.e.:

a = mt^2+nt+r

b = pt^2+qt+s

Really thankful for all your help!

[Torsten]

What's the problem ?

function dxydt = rate(t,xy);
  m = ...;
  n = ...;
  r = ...;
  p = ...;
  q = ...;
  s = ...;
  a = m*t^2+n*t+r;
  b = p*t^2+q*t+s;
  x = xy(1);
  y = xy(2);
  dxydt = [a-b+c*x*y;c*y*(1-y)-a*y/x];
end