# Solving a system of differential equations

I have a system of 5 differential equations with 5 unknown variables So I have 4 equations differentiated with respect to time and the 5th equation is a partial differential equation with respect to time and distance. All the equations have variables interdependent with each other. Example:

   dF/dt = 3*F  +  4*G  +  Constant;
dG/dt = 3*H  +  4*G  +  Constant + F;
dH/dt = 3*H  +  4*G  +  Constant  +  I;
dI/dt = 3*H  +  4*I  +  J*5  +  Constant;
∂J/∂t = 3*F  +  4*H  +  J  +  ∂J/∂x;


I am solving this in matlab. I have initial values and boundary conditions and it is differentiated for a period of 24 hours .I tried solving using explicit method,encountered a lot of iteration errors I have also tried solving using implicit finite difference method.But this is not giving me any good results. I have also tried using differential syntax to solve them. I heard about time series,and how does it apply here if possible?Are there any other methods I can approach this system?

• Is there anything important about the last line not being $d/dx$? – Simply Beautiful Art Aug 27 '16 at 20:06
• the d/dx is essentially fluid flow(motion) that is in the x direction with respect to time while the rest of the equation variables are static – John Aug 27 '16 at 20:19
• Two things : 1) If $\partial J/\partial x$ having small evolutions compared to $\partial J/\partial t$. In this case, you could consider it as constant, solve your system for a small interval of time, reajust $\partial J/\partial x$, etc. 2) You say you have used "differential syntax": doyou mean by that ode45 for example ? Much better than re-implementing explicit or implicit methods But there is better to do, it is by using expm (matrix exponential) ; of course, constants have to be incorporated... – Jean Marie Aug 27 '16 at 20:39
• Thank you for the suggestion!greatly appreciated. I did try the first method as increments in x are very small compared to t.Exponential method looks promising. However if the constant method is not considered then how does solving by ode work?Is it possible to solve the first four using ode and the last one by pedpe? – John Aug 27 '16 at 20:56
• @JeanMarie I have tried to use the exponential method but i am not sure how to go about it – John Aug 29 '16 at 9:28

This is not a solution. Just a slightly elaborate approach.

Let us consider, for the explanation of the idea, an (over)simplification :

• all constants are $0$, and as well, (see our discussion)

• $\partial J/\partial x=0$ (considered as a constant).

If $A:=\pmatrix{3&4&&&\\1&4&3&&\\&4&3&1&\\&&3&4&5\\3&&4&&1}$, the solution is $\pmatrix{F(t)\\G(t)\\H(t)\\I(t)\\J(t)}=exp(tA).\pmatrix{F(0)\\G(0)\\H(0)\\I(0)\\J(0)}$

In Matlab (I rely on my memory), it would be implemented more or less like that (using "expm" function, not "exp"):

A=[3 4 ...; ... ; 3 0 4 0 1];

V0=[1,3,0,0,0]'; % or any other set of initial values.

n=1000;

T=zeros(5,n);

for k=1:n

t=k/100; % time

V=expm(t*A)*V0;% V0 has to be a column vector

T(1:5,k)=V;

end

plot((1:n)/100,T(1,:)) ; % example : the evolution of variable F

Important remark: The eigenvalues of $A$ are

$$\{8.2, 3.69 \pm 0.8i, -0.29 \pm 1.48i\}$$

It attracts attention on the fact that there is a trend to "global" exponential growth (the dominant eigenvalue $8.2$ is $>0$) with periodic fluctuations.

Edit

One can improve the computational efficiency without impairing the overall precision by using the fact that :

$$\tag{1}\forall t1,t2, \ \text{expm}(t1*A)*\text{expm}(t2*A)=\text{expm}((t1+t2)A)$$

(Remark: it is not true in general that expm$(M_1)*\text{expm}(M_2)=\text{expm}(M_1+M_2)$; but it is true in particular when $M_1$ and $M_2$ can be expressed as polynomials in the same matrix, as is the case here).

If relationship (1) is applied with $t_2=t$ and $t_1=\Delta t$, this gives rise to the fact that it suffices to left-multiply the current state $\text{expm}(t*A)*V_0$ by a multiplicative factor, always the same : $M=\text{expm}((\Delta t)*A)$ to obtain the new current state.

This gives rise to the following modification of the heart of the previous program:

V = V0;

Deltat = 0.01;

M=expm(Deltat*A);

for k=1:n

T(1:5,k)=V;

V=M*V;

end

• Have you seen my answer ? – Jean Marie Aug 29 '16 at 13:31
• Yes, thank you greatly appreciated , i was working on it for 5 hours, just got it to work from your method, it took some time as my equations were a little more complicated. I have voted your solution as answer, Cheers. – John Aug 29 '16 at 17:51
• Happy that it has helped you a little. The program I have given (in which I just corrected a little flaw) can be greatly improved. If you are interested, I can indicate you how. I imagine you are working on some aspects of chemical cinetics ? – Jean Marie Aug 29 '16 at 17:57
• That would be brilliant.I am working on modelling by energy balance equations. – John Aug 29 '16 at 20:06
• I will add something as an edit to my answer. – Jean Marie Aug 29 '16 at 20:28