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
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.
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
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
on 15 Mar 2023
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
on 15 Mar 2023
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
on 16 Mar 2023
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
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
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
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
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!
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
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
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
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
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
Find more on Numerical Integration and Differential Equations in Help Center and File Exchange
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 (한국어)