**You are now following this question**

- You will see updates in your followed content feed.
- You may receive emails, depending on your communication preferences.

# Please help solve a system of differential equation

4 views (last 30 days)

Show older comments

### Answers (1)

Alan Stevens
on 15 Mar 2023

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

##### 19 Comments

Jayesh Jawandhia
on 15 Mar 2023

Edited: Jayesh Jawandhia
on 15 Mar 2023

Alan Stevens
on 15 Mar 2023

Jayesh Jawandhia
on 15 Mar 2023

Alan Stevens
on 16 Mar 2023

Alan Stevens
on 16 Mar 2023

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
on 16 Mar 2023

Alan Stevens
on 17 Mar 2023

- What is the system you are trying to calculate? i.e. What do x and y physically represent?
- How did you derive the equations, and why does your data set contain estimates for a and b, but not c?
- 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
on 17 Mar 2023

Hey Alan,

Trying to answer your questions below:

- 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.
- We derived these equations using mass balance (Inlet - Outlet + Generation = Accumulation) on both air and liquid.
- 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
on 18 Mar 2023

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

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

Jayesh Jawandhia
on 18 Mar 2023

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
on 18 Mar 2023

Sam Chak
on 19 Mar 2023

Edited: Sam Chak
on 19 Mar 2023

@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
on 19 Mar 2023

Edited: Jayesh Jawandhia
on 19 Mar 2023

Sam Chak
on 19 Mar 2023

Edited: Sam Chak
on 19 Mar 2023

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
on 19 Mar 2023

@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
on 19 Mar 2023

Jayesh Jawandhia
on 4 Apr 2023

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
on 4 Apr 2023

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

### See Also

### Categories

### Community Treasure Hunt

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

Start Hunting!**An Error Occurred**

Unable to complete the action because of changes made to the page. Reload the page to see its updated state.

Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list

How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

Americas

- América Latina (Español)
- Canada (English)
- United States (English)

Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)

Asia Pacific

- Australia (English)
- India (English)
- New Zealand (English)
- 中国
- 日本Japanese (日本語)
- 한국Korean (한국어)