Solving Coupled Differential Equations

For my problem I have two differential equations. These are:

dx/dt = -xp + yq + (1-x-y)q
dy/dt = -yq + (1-x-y)u

For my problem I have a given set of historical data and my objective is to find the best values for p,q & u to give the best fit to this historical data. I have a column of historical data related to each equation above. I also need to find the best initial conditions for each equation.

So far I have attempted to use For loops for each of p,q&u in order to find the best values compared with my data.

My problem is that I am not too sure on how to code this correctly. If anybody could help me make a start that would be great. I understand the general idea but am unsure of how to carry this out in Matlab

[Shyamal Laxman Hanabar-]

Write the matlab code to solve integration of (3*x+4)/(4x+5)dx

[Shyamal Laxman Hanabar-]

Write matlab code for find du/dz if u=loge(x^3+y^3+z^3-3*x*y*z^3)

[Rik]

I moved your (@Shyamal Laxman Hanabar-) comments to the question, since they are not answers. You have already received advice on your separately posted question.

Answer 1

[Star Strider]

If you are fitting a system of differential equations to data, see Monod kinetics and curve fitting. It uses estimated parameters to define the initial conditions as well as fit the data. It should be straightforward to adapt it to your problem.

[Dave Black]

Thank you for getting in touch over this. I have had a look at the link but am still a little unclear. Perhaps I could provide more information about my problem? For example I do know the range that my p,q & u values fall within. And I do know the range for which my initial x & y fall within. But I am having issues with determining a good solution compared to my known data. It would be great if you could provide anymore guidance.

[Star Strider]

My pleasure.

The initial parameter estimates are important, so you may need to experiment with several different sets before you arrive at the correct set. Also, the Monod kinetics problem only output as a vector in ‘S’. You apparently need to output and fit both ‘x’ and ‘y’, so you need to fit a (Nx2) matrix of ‘[x(:) y(:)]’ values that you need to create and then fit, so you pass the (Nx2) matrix as a single dependent variable to lsqcurvefit. The lsqcurvefit function can do that, and you have to write your objective function to accommodate it. Also, your ‘tspan’ needs to be the time vector of your data.

[Dave Black]

Thanks a lot for responding again. I would like to go into a little more detail about this with yourself. Basically the idea that your given approach follows is different to what I have had suggested to me. I have been advised to follow loops for all possible values of p q & u and x & y initial conditions firstly. Doing this in a loop will give me many many answers as there are so many different possibilities. This is over a timescale of 0-6. I have moved over my excel data where I have two columns and of course solving my equations will give me two lines when I plot. I have been advised to somehow start with a large sum figure. Then work out the absolute difference between the plots given and my known data points. Initially of course this value will be larger than my sum figure I have provided. So effectively after running one loop the answer of plot-data points will become the new sum. The idea is to keep rerunning the loop until I have the smallest difference between my equation solution and data points. Would you have any knowledge on how to implement this? I know it is difficult to get to grips with over plain text but I would be happy to provide more information.

Thanks again for responding

[Star Strider]

‘I have been advised to follow loops for all possible values of p q & u and x & y initial conditions firstly.’

That advise is simply wrong. The entire Optimization Toolbox, and all the theoretical work that preceded it, is designed to make parameter estimation both intelligent and efficient. The method you have been advised to use will take years to converge, if it ever converges on the optimal values, which would be sheer luck.

Using the Optimization Toolbox functions such as lsqcurvefit will probably take a few minutes at most to converge on the best-fitting parameter estimates. (The Global Optimization Toolbox has other functions that will search more widely for the best parameter estimates, if lsqcurvefit has problems.)

You can write your own genetic algorithm if you do not have the Optimization or Global Optimization Toolboxes. It will take longer than the optimised functions in the Toolboxes, but it is infinitely more efficient than the ‘random-luck’ search you have been advised to use.

I will help you parameterise your differential equation ODE function if you need help with it. (Have you already tested it to be sure that it works with ode15s?) I do not have your data, so I cannot help you with the code to fit it to your ODE system beyond what I have already advised.

I will help further if I have more information. I am reluctant to create my own functions and data because that is rarely efficient for such problems.

[Dave Black]

Thank you for responding again, is there a way I can contact you in order to give you more information such as data values?

I have a system working where I know the values for p q & u and when I provide them it works to produce a graph that fits my data. I have used ODE23 to solve.

I can certainly provide more detail and show my work so far if you could help somehow?

Thanks!

[Star Strider]

My pleasure.

I prefer to keep everything here, for continuity. That makes this thread (and other such threads) valuable for others with similar problems.

The fitting procedure is not difficult. I need your data (or a representative sample of it, preferably as a ‘.mat’ file) to write code specific to it for lsqcurvefit. It would help significantly to have the differential equation function you wrote so I can parameterise it in context.

[Dave Black]

That would be really helpful. Is there a way I can pass on a .mat file to yourself on this thread? I am not too familiar with this site. So far I have stored my known data as a matrix on matlab and have attempted to use ode23 to solve my equations but I can show you this no problem.

I am in the U.K. so I can send the file in the morning if that would suit? Do not worry if you are outside the UK I am in no real rush for immediate help. But anything you could do would be really beneficial.

Thanks!

[Star Strider]

My pleasure.

Create a ‘.mat’ file with the save function and save your ‘time’, ‘x’, and ‘y’ data to it. Then use the ‘paperclip’ icon (next to the green and brown ‘picture’ icon) to upload it.

You can copy the code of your ODE function and paste it to a Comment here. Please highlight your code, then use the {}Code button to format it.

It would be convenient if you copy and paste your code to the same Comment you attach your data to.

[Dave Black]

That's great I shall do that in the morning. Hopefully all that I send over will make sense to yourself. I think the problem sounds more complex than it actually is.

Thanks again!

[Star Strider]

My pleasure!

[Dave Black]

Hello Star Strider and thanks again for offering your help. I have put together my sample data and attached it here in the form of a .mat file. Hopefully you can open it and the script file that I have also attached.

Here is the code I have been attempting to use so far, I know this code works when I know what p,q & u are as well as knowing my initial conditions.

I firstly created a function file like so:

function [Ddv_Div] = ReligionAffiliation(I, D)

p = 0.070;
q = 0;
u = 0.05;

x = D(1);
y = D(2);

Ddv_Div = [ -x*p + y*q + (1-x-y)*q;
    -y*q + (1-x-y)*u ];

end

Then I ran this function within a script file :

domain = [0 6];

IC1 = 0.47;
IC2 = 0.03;

IC = [IC1 IC2];

[IVsol, DVsol] = ode23('ReligionAffiliation' , domain, IC);

plot(IVsol, DVsol(:,1), 'k') % x
hold on
plot(IVsol, DVsol(:,2), 'r') % y

As you can see from this code I have used my known figures for p,q,u & initial conditions. The "domain" of 0 to 6 is essentially the same as giving my years "1951 1961 1971...".

A few points to note: For my problem I have 6 different models that consider religious affiliation, the code for the differential equations I have given is for Model 1 which is all I am worried about solving for now. I would like to know what the best values are for p q & u to fit the data I have provided as well as the best initial starting conditions.

However, for all 6 models I do know the ranges that each variable will fall within. For example I know that p will always be between 0 & 0.6, q will be between 0 & 0.01 and u will always be between 0 & 0.5.For initial X starting I know it'll be between 0.4 & 0.5 and for Y it will be between 0 & 0.05.

My data does not change based on the model.

I have also been advised to forget about the "1981" data point when attempting to find the best fit because it is an anomaly.

I hope I have provided enough information and as mentioned any help at all would be greatly appreciated.

[Star Strider]

That should be enough. It will take some time, so I will work on this and post back here later. (I will also delete this Comment.)

[Dave Black]

That would be perfect! I do know the best answers for model one as given in the research paper I am following. So p = 0.070,q = 0,u = 0.050 and the initial conditions are around x = 0.47 and y = 0 to 0.02. So these figures can be used to know if I am on the right lines

Thanks again!

[Star Strider]

My pleasure!

I created a function file (with no arguments or outputs), put all my code here in it, and then just ran the function. It is much easier than running everything from a script.

The Code —

D = load('ReligionData.mat');
RD = D.religiondata;
t = RD(:,1);
xy = RD(:,2:3);

      function X = odefit(b,t)
      % % % b(1:3) = Parameters,  b(4:5) = Initial Conditions
      % % % MAPPING: x(1) = x,  x(2) = y,  b(1) = p,  b(2) = q,  b(3) = u
      Ddv_Div = @(t,x,b) [ -x(1).*b(1) + x(2).*b(2) + (1-x(1)-x(2)).*b(2);  -x(2).*b(2) + (1-x(1)-x(2)).*b(3) ];
      [~,X] = ode45(@(t,x) Ddv_Div(t,x,b(1:3)), t, b(4:5));
      end

X0 = [0.07;  0;  0.05;  1;  1];

B = lsqcurvefit(@odefit, X0, t, xy, zeros(size(X0)), inf(size(X0)))

fprintf(1, '\n\tParameters:\n\t\tp = %11.3E\n\t\tq = %11.3E\n\t\tu = %11.3E\n\tInitial Conditions:\n\t\tx = %11.3E\n\t\ty = %11.3E\n', B);

Ddv_Div = @(t,x,b) [ -x(1).*b(1) + x(2).*b(2) + (1-x(1)-x(2)).*b(2);  -x(2).*b(2) + (1-x(1)-x(2)).*b(3) ];
[T,X] = ode45(@(t,x) Ddv_Div(t,x,B(1:3)), t, B(4:5));

figure(1)
plot(t, xy, 'p')
hold on
plot(T, X, '--')
hold off
grid
legend('Data: x', 'Data: y', 'Fit: x', 'Fit: y', 'Location','W')

Producing:

  Parameters:
    p =   2.817E-02
    q =   6.515E-03
    u =   1.370E-06
  Initial Conditions:
    x =   9.217E-12
    y =   4.256E-01

I initially got negative initial conditions, so I put constraints on the parameters to prevent that. With the constraints, the initial condition for ‘x’ is essentially 0, and for ‘y’, 0.43.

The Plot —

[Dave Black]

Thanks so much for this this definitely seems along the right lines towards what I am after. Just a couple of questions I have, I'm just wondering where you have X0 = [0.07;0;0.05;1;1] I believe these are the p q and u values I mentioned previously. I'm just wondering how this is used because these are the actual answers I want for p q and u where in the end you have p=0.02,q=0.006,u=0.000001

I understand it's a little difficult when you don't have a greater understanding of the topic but basically I am following a specific research paper so I know the answers that provide the best fit are p=0.07,q=0 and u =0.05 . But I am supposed to produce a system that finds these answers as if I don't know them.

I hope this is clear! The fits do of course look correct graphically so I am very grateful for this method and it is a great help.

[Star Strider]

My pleasure.

I used your values as initial parameter estimates (the last two were initial condition estimates), then let lsqcurvefit find the parameters that provided the best fit. I do not know the parameter estimation methods the authors of the paper you alluded to used.

I am confident that the values here are accurate. I checked the code to be certain the values returned by the objective function matched the values in your data (so ‘x’, or ‘x(1)’ is in the first column, and ‘y’ or ‘x(2)’ is in the second). I trust the parameter estimates lsqcurvefit returned. The graphical fit is quite good

I would not be completely trusting of the information in published articles, even in refereed journals. That a paper is published does not mean it should have been. I have caught several reported models in my career that had erroneous parameter estimates, since either the authors used some sort of unjustified linearization technique (doing a linear fit of the logarithm of an exponential function is a frequent error committed by those who do not understand or have access to nonlinear parameter estimation techniques, as is logarithmic ‘curve-stripping’ to estimate the parameters of a multiple-exponential model, these — and others — producing wildly incorrect parameter estimates), or the authors simply did not understand regression techniques. Some were so clearly wrong that I wonder how such obvious errors got past the reviewers. One paper attempting to model glucoregulation in normal subjects estimated 15 parameters from 13 data pairs, and could not explain the reason the confidence intervals on the parameters all included 0! (I have been a reviewer, and caught several such errors in manuscripts, and have helped authors correct them.)

Being unable to reproduce published results is not at all uncommon. I trust my results.

If my Answer helped solve your problem, please Accept it!

[Dave Black]

I understand now why you have supplied those values, they are sort of an initial guess as to what you think the values will be?

As far as I know the author of the paper used a Brute force approach and it was all coded on fortran, I do have all of this code but can't follow it on Matlab unfortunately.

You clearly seem to have a strong understanding of what you are talking about so I am very grateful for your help and it will be very beneficial to myself.

I only have one more query,say for example I wanted to forecast the plots into the future so for example stretch the line plots past the data into the year 2151 is there a way I could do this?

Again I do appreciate the time you have taken to help!

[Star Strider]

The values lsqcurvefit returned are the parameter estimates it arrived at using its optimization algorithms. The initial parameter estimates I used to give lsqcurvefit necessary initial parameter estimates were the ones you provided. The ones lsqcurvefit estimated were quite different. I also did not find it necessary to eliminate any data, since lsqcurvefit and my code were able to use all of them.

A ‘brute force’ approach in FORTRAN or any other language is no substitute for efficient optimisation algorithms that have been available in some form for at least the last 50 years. I coded a Lavenberg-Marquardt nonlinear parameter estimation algorithm in Apple FORTRAN in 1982. A ‘brute force’ approach is never appropriate. I would not trust the parameter estimates of a ‘brute force’ approach.

Please do not attempt to extrapolate that far beyond the region of fit. It is theoretically possible to ‘get numbers’, but they would be meaningless. You would first need to calculate confidence intervals on the fit (not possible with the Statistics and Machine Learning Toolbox functions with lsqcurvefit using two dependent variables), so you would have an estimate of the prediction error. You could see if you could integrate them analytically and fit each set of data separately using the parameters my code estimates with nlinfit and then use nlpredci to calculate the confidence intervals. I suspect that the confidence intervals on the prediction would expand to include 0 long before you reached 2151, making the estimate itself worthless. (For fun, the values I get for 2151 are 0.1877 and 0.1158 for ‘x’ and ‘y’ respectively. The curves cross at about 2075.)

As always, my pleasure.

[Dave Black]

I really do appreciate your wisdom on all of this and I will attempt to apply the methods you have put forward.

Thank you for all you have done.

[Star Strider]

My pleasure.

You can use your other models with my code. If they are in ‘x’ and ‘y’ only, parameterise them as I did with this model, and you can use them with only the necessary changes in the parameter vector length and parameter mapping.

[Dave Black]

I hope you get to see this comment. I have tried to implement your code on Matlab but I think that I am getting the ordering wrong somewhere and it is not wanting to run for me, also how did you apply the non negative contsraints to this?

[Star Strider]

If you copy and paste my entire posted code to a function file (as I did), it should run without error. (I’m running it in R2017a, but I doubt there are any differences in recent versions.)

If you’re modeling different differential equation systems than the one in my code, the order in which you enter them into the ODE function (here ‘Ddv_Div’) is extremely important. The first row will be returned as the first column in the ode45 solution matrix, and the second row in the second column. The data you’re fitting have to be columns in your dependent variable matrix, and match those column assignments, or the lsqcurvefit function will return erroneous results for the parameter estimates.

I set the constraints in the lsqcurvefit call here:

B = lsqcurvefit(@odefit, X0, t, xy, zeros(size(X0)), inf(size(X0)))

setting the lower limit, ‘lb’ in the documentation, as ‘zeros(size(X0))’, and ‘ub’ as ‘inf(size(X0))’.

If negative parameter values are permitted, omit both ‘lb’ and ‘ub’ from the ‘lsqcurvefit’ arguments. If some are constrained and others aren’t, set up the specific elements of those arguments to the limits you want. Note that -Inf and Inf are both acceptable constraints. See the documentation for the lsqcurvefit function for details.

[Dave Black]

Hello Star Strider,

Again I hope you get to see this comment.

Your method has been the best I have come across for solving this problem, however, I would still like to extrapolate the lines out further into the future for some general analysis. Could you provide any advice on how to code this sort of thing?

[Star Strider]

Thank you!

Extrapolating is straightforward. Just extend the upper limit of the time vector.

To extrapolate to 2500, this works:

tspan = [min(t) 2500];
Ddv_Div = @(t,x,b) [ -x(1).*b(1) + x(2).*b(2) + (1-x(1)-x(2)).*b(2);  -x(2).*b(2) + (1-x(1)-x(2)).*b(3) ];
[T,X] = ode45(@(t,x) Ddv_Div(t,x,B(1:3)), tspan, B(4:5));

figure(1)
plot(t, xy, 'p')
hold on
plot(T, X, '--')
hold off
grid
legend('Data: x', 'Data: y', 'Fit: x', 'Fit: y', 'Location','W')

I do not suggest extrapolating beyond the region of fit, but it can be an interesting exploration.

[Dave Black]

Hi Star Strider,

Thanks for this. Realistically I should never have to extrapolate beyond 2250 so this should be fine. It is just to note when the top line (X related line) drops to 0.1 and the year when Y (bottom line) gets to 0.5. I had to switch the data columns around for the results to make sense in the real world. It's strange though, my advisor who used brute force has results and mine are all very similar but out by a factor of 10. At least it all finally makes sense though.

[Star Strider]

My pleasure.

If you and your advisor are using the same code for the differential equation, and the same parameters, I cannot explain the results being off by a factor of 10. The data columns must match the ‘rows’ of ‘Ddv_Div’, so the first column matches the first row, the second column the second row.

Personally, I trust the results of the optimised lsqcurvefit parameter estimation here. I would not put much stock in any brute force method.

[Dave Black]

It's just odd isn't it. Yes, I switched so they match and now the results do indeed make sense. Also when I extrapolate the graphs out the analysis results I am looking for (e.g. the year when x<0.1) match perfectly those of my adviser. So I think this is it working and I thank you so much for this

[Star Strider]

My pleasure.

Using this optimised parameter estimation, your results are the correct results!

Answer 2

[Roger Stafford]

I would think it would be far more efficient to do your fitting using a general solution to your two linear differential equations. That way you could use some of the fitting tools that matlab has. Happily, according to dsolve (after a little simplification) the general solution involves just exponential functions and is:

x = C1*(p-u)/u*exp(-(p+q)*t) + q/(p+q);
y = C1*exp(-(p+q)*t) + C2*exp(-(u+q)*t) + p*u/((u+q)*(p+q));

where p, q, and u are the given parameters and C1 and C2 are dependent on initial values of x and y. This means you have a straight-forward five-dimensional fitting problem to solve, rather than having to deal with differential equations.

Note: Are you sure about the first differential equation

dx/dt = -x*p + y*q + (1-x-y)*q  ?

Notice that the y terms cancel each other, giving the simpler

dx/dt = -(p+q)*x + q

expression which does not involve y.