Solving a system of differential equation numerically but with constants

Hello,

I am asking this question before I buy MATLAB to make sure that I can get a solution of a problem (so that I won't be wasting money). I have never used MATLAB before.

I am currently needing a numerical solution to a system of differential equation for a certain phenomenon I am currently working on. However, the problem is when I look through google, most question involves a system of differential equation with constants that are defined apriori (like 3*x^2 instead of a*x^2). I am using the system of differential equation to fit certain experimental data I have achieved so that I can figure out the constants that have important meaning to the phenomenon. Therefore, I have no idea what value the constants take.

The system of differential equation I have set up looks like this:

w(t) = a*exp(-t/b)

dx(t)/dt = c*w(t) + d*y(t) + e*z(t) - (1/f)*x(t)

dy(t)/dt = g*x(t) + h*y(t) + i*z(t) - (1/j)*y(t)

dz(t)/dt = k*x(t) + l*y(t) - (1/m)*z(t)

where x(0) = y(0) = z(0) = 0, dx(0)/dt = dy(0)/dt = dz(0)/dt = 0, and all the constants from a to m is above 0.

There are four differential equation and four yet-to-be-known equation, so I think this is theoretically solvable. However, each equations are recursive (am I using this term right?) and I cannot do this by hand and nor could a freeware software that tried to solve this analytically.

I am expecting a solution where each equation consists of variable (t) and constants combined in an arithmetic way. Is this possible in MATLAB?

I really have no background in mathematics. Even worse, I learned almost everything in Japanese so I do not know the correct translation of tehcnical terms. Also, my English is mediocre at most. I apologize for my inability to effectively explain what I am trying to explain.

Thank you.

[Torsten]

In principle it is possible with MATLAB.

You will have to couple a solver for ordinary differential equations (e.g. ODE15s) with a nonlinear least-squares solver (e.g. LSQCURVEFIT from the optimization toolbox).

The nonlinear least-squares solver suggests values for your constants, the ODE solver returns the solution at the times where you have data from your experiments for these constants, the nonlinear least-squares solver tries to adjust the constants such that x, y and z better meet your experimental data and so on in an iterative fashion.

Take a look at

 http://www.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting

where a similar problem has been discussed.

But parameter fitting problems are inherently difficult. Nobody can guarantee that MATLAB (or any other code) will produce good results. Much depends on your experimental data and whether your model is adequate or not.

I suggest you try a MATLAB test version.

I heard that TOMLAB is especially suited for parameter fitting problems, but I myself don't have experience with it.

What makes me wonder in the formulation of your mathematical problem is that you prescribe six boundary conditions for three ordinary differential equations. This does not seem to make sense to me. Think about your problem formulation before you start with numerical computations.

Best wishes

Torsten.

[Torsten]

By the way:

If you write the equation for w as

dw/dt = -1/b*w(t),

you get a linear system of 4 ODEs with constant coefficients which can be solved analytically.

Best wishes

Torsten.

[Hayao]

Thank you Torsten!

Sounds like it's possible. Great! I'll probably have tough time trying to get used to the interface and be able to effectively use it though...

About the boundary conditions, some of those boundary conditions were not necessary. I just listed 6 even though I just needed 3 because I already know that it is going to be like that one way or another (due to the property of the phenomenon).

Also, you commented about how they can be analytically solved if I write the equation for w as you mentioned. I tried it on Maxima, a freeware computer algebra system, and it came out with inverse Laplace transformed equations that are so complicated to a degree where they are incomprehensible.

Am I doing something wrong? Or is it a software problem?

[Torsten]

The solution is "simply" exp(A*t) where A is the coefficient matrix:

http://mathworld.wolfram.com/OrdinaryDifferentialEquationSystemwithConstantCoefficients.html

Best wishes

Torsten.

[Hayao]

Okay thank you for the advice. I decided to do this by hand as much as possible. Correct me if I am wrong because this is the first time I am doing something like this. I am supposed to get a solution as a system of equations expressed in vector:

x(t) = C exp( A t)

where C is the value for each equation when t=0 expressed in vector, and A is a matrix of constants used in each differential equation. I am supposed to solve exp( A ) so that I can achieve each coefficients of the solution equations. If I am not mistaken, this can be solved as matrix exponential:

exp( A ) = SIGMA(k = 0 to infinite) [( A ^k)/k!]

This should theoretically converge since it's a 4x4 square matrix. But I cannot do this by hand in the case I am working on because the matrix A is defined as:

A =

[ a11 0 0 0 ]

[a21 a22 a23 a24]

[a31 a32 a33 a34]

[ 0 a42 a43 a44 ]

This is not solvable by hand. I tried using Maxima to compute the matrix exponential in hope of achieving a converged solution, but it was unable to.

Can this be solved in MATLAB?

[Torsten]

The matrix exponential of a (4x4) general matrix will be a lengthy expression and it will involve the roots of a polynomial of degree 4.

Take a look at

http://www.wolframalpha.com/input/?i=matrix+exponential&a=*C.matrix+exponential-_*Calculator.dflt-&f2=%7B%7Ba%2C+0%2C+0%2C+0%7D%2C+%7Bb%2C+c%2C+d%2C+e%7D%2C+%7Bf%2C+g%2C+h%2C+k%7D%2C%7B0%2Cm%2Cn%2Co%7D%7D&f=MatrixExponentialCalculator.matrixexp_%7B%7Ba%2C+0%2C+0%2C+0%7D%2C+%7Bb%2C+c%2C+d%2C+e%7D%2C+%7Bf%2C+g%2C+h%2C+k%7D%2C%7B0%2Cm%2Cn%2Co%7D%7D

for a solution.

Best wishes

Torsten.

[Hayao]

Okay, so that's why it's not easily solved.

That's good to know. I'll try more simplification to my model so that the exponential matrix can be solved more easily.

Thank you Torsten for very helpful comments. I never imagined that it is possible to solve a system of differential equation that is relatively complicated to solve by hand using matrices. New knowledge!

[Torsten]

I think there is no need to simplify your model.

Solving your problem numerically, your constants will always be definite numbers. So a simple call to MATLAB's matrix exponential as expm(A*t) will solve your ODE system.

The solution in MATHEMATICA only looks so complicated because the entries of the matrix are given as arbitrary parameters.

Best wishes

Torsten.

[Hayao]

The constants are arbitrary. I do not know what number they actually are. I am not exactly trying to solve the ODE system. Instead, I am trying to find out what the constants are by fitting the ODE system to an experimental result.

As a matter of fact, it is virtually impossible to expect what the constants are. Some of them can be approximated, but lack of experimental method keeps those constants to be unknown (there is a way, but I don't have the machine capable of doing it).

[Torsten]

Yes, of course a,b,c,.. are unknowns, and they will change during the computation until the best fit for your model is reached. But if you solve your problem with MATLAB, a,b,c,... will always be (changing) numerical values such that you apply the matrix exponential to a matrix with definite numbers as entries, not to a symbolic matrix with letters as entries.

Best wishes

Torsten.

Answer 1

[Gajendra Singh]

Do you find any solution?? Actually i have same problem ( four differential equation with 4 unknown variable and 4 unknown constant which have to be optimize)