Hello community
I need to solve a complex set of equations using Mathematica and would like to know if you could help with the writing in the program. These are the equations that need to be written:
dx/dt=-ax-bxy-dxz
dy/dt=ax-by+cz-dxy
dz/dt=by-cz-cxz
d?/dt=dxy+exz
where x,y,z,? are variables and a, b, c, d, e are constants.
How can I obtain an expression for each variable that depends on the constants? Which are the functions to introduce in Mathematica?
Thanks
Your first step is to learn about Mathematica. Take a look at some of the references listed in the following link and start to learn the system, then take a try at your problem. The particular function that you will be interested in is DSolve.
DSolve did not give solution:
ClearAll[x, y, z, gamma, t, a, b, c, d, e];
eq1 = x'[t] == -a x[t] - b x[t] y[t] - d x[t] z[t];
eq2 = y'[t] == a x[t] - b y[t] + c z[t] - d x[t] y[t];
eq3 = z'[t] == b y[t] - c z[t] - c x[t] z[t];
eq4 = gamma'[t] == d x[t] y[t] + e x[t] z[t];
DSolve[{eq1, eq2, eq3, eq4}, {x[t], y[t], z[t], gamma[t]}, t]
But NDsolve can solve them. You need to supply numerical values for all the parameters and initial conditions also.
ClearAll[x, y, z, gamma, t, a, b, c, d, e];
a = 1; b = 2; c = 3; d = 4; e = 5;
eq1 = x'[t] == -a x[t] - b x[t] y[t] - d x[t] z[t];
eq2 = y'[t] == a x[t] - b y[t] + c z[t] - d x[t] y[t];
eq3 = z'[t] == b y[t] - c z[t] - c x[t] z[t];
eq4 = gamma'[t] == d x[t] y[t] + e x[t] z[t];
NDSolve[{eq1, eq2, eq3, eq4, x[0] == 1, y[0] == 0, z[0] == 0,
gamma[0] == 0}, {x[t], y[t], z[t], gamma[t]}, {t, 0, 10}]
Thank you Nasser. Is there any way to obtain a formula corresponding to this interpolated functions or plots?
Is there any way to obtain a formula corresponding to this
interpolated functions or plots?
I am no expert in this area, but the way I’ve seen it done is by getting the x,y values from the InterpolatingFunction object, then use InterpolatingPolynomial.
InterpolatingFunction itself uses piecewise InterpolatingPolynomial under the cover. But to obtain just one InterpolatingPolynomial over all the grid points, will result in a very large polynomial of course (order one less that number of grid points), and is not recommended as it will be highly oscillatory. You could pick smaller grid regions at a time, and find the InterpolatingPolynomial over those. But here is what you can try:
ClearAll[x, y, z, gamma, t, a, b, c, d, e];
a = 1; b = 2; c = 3; d = 4; e = 5;
eq1 = x'[t] == -a x[t] - b x[t] y[t] - d x[t] z[t];
eq2 = y'[t] == a x[t] - b y[t] + c z[t] - d x[t] y[t];
eq3 = z'[t] == b y[t] - c z[t] - c x[t] z[t];
eq4 = gamma'[t] == d x[t] y[t] + e x[t] z[t];
sol = First@NDSolve[{eq1, eq2, eq3, eq4, x[0] == 1, y[0] == 0, z[0] == 0, gamma[0] == 0},
{x, y, z, gamma}, {t, 0, 10}];
xSol = x /. sol; (* get formula for the x(t) solution *)
xpoints = Flatten@xSol["Coordinates"];
ypoints = Flatten@xSol["ValuesOnGrid"];
data = Transpose[{xpoints, ypoints}];
p = InterpolatingPolynomial[data, t];
Expand[p]
Here is your formula for the x(t) solution found by NDSolve
