Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm trying to fit my data to a system of ODEs:

x'[t]==a*x[t]+b*y[t]
y'[t]==c*y[t]+d*x[t]+e*z[t]
z'[t]==f*z[t]+g*y[t]

having t,x[t],y[t],z[t]. I followed the code from here, but in the end got an error "Part::partw: "Part 2 of {{X[0]->6.99188*10^7,Y[0]->1.44163*10^6,Z[0]->1.22418*10^6}} does not exist. General::stop: Further output of Part::partw will be suppressed during this calculation. >>"

I'm a newbie in mathematica, can you point out what I do wrong? Or maybe you can suggest better way to do such fitting? Here's my code:

data = {{0.0001, 69721851.824697, 1710041.15718086, 1249648.07662318},
{0.0002, 69043570.7566936, 2066441.5110575, 1844749.71071548},
{0.0005, 67932147.9362317, 2994541.62071368, 3108260.41352905},
{0.0009, 66863760.1248742, 3587039.29801347, 4282635.21874863},
{0.0015, 64345357.8511309, 968094.021795217, 10585640.3555617},
{0.003, 64768407.2504181, 3655819.53710869, 7812897.34991382},
{0.005, 63761257.6277216, 1222706.2166919, 12090416.0146563},
{0.01, 64307991.5530707, 1139370.28069903, 12657324.2458834},
{0.05, 68583353.8471158, 1802661.70625935, 13079855.4358155},
{0.08, 71169929.0580033, 1496318.92452295, 13991708.7133703},
{0.1, 72811456.5338936, 1853138.14008422, 13681128.3978895},
{0.2, 78742983.529449, 3190647.61480991, 13511961.6223572},
{0.5, 88513835.0051046, 3556145.38852813, 14771090.342085},
{1, 94765316.2395368, 4882322.77031356, 14796960.6588559},
{2, 98437836.0564595, 5353800.46123903, 14862252.8596808},
{5, 99708229.0609698, 3726362.20580859, 16993877.2198709},
{10, 99634510.8439181, 4858332.59695385, 16478819.0906012}};
{t, x, y, z} = 
 Transpose[data];(*extract data*){xdata, ydata, zdata} = (Transpose[{t, #}]) & /@ {x, y, z};(*get {t,x} and {t,y} pairs*)fx =  Interpolation[xdata, Method -> "Spline"];(*spline interpolation*)fy =  Interpolation[ydata, Method -> "Spline"];
fz = Interpolation[zdata, Method -> "Spline"];
dx = fx'[t];
dy = fy'[t];
dz = fz'[t];
ab = FindFit[Transpose[{fx[t], fy[t], dx}], a*X + b*Y, {a, b}, {X, Y}]
cde = FindFit[Transpose[{fx[t], fy[t], fz[t], dy}], c*Y + d*X + e*Z, {c, d, e}, {X, Y, Z}]
fg = FindFit[Transpose[{fz[t], fy[t], dz}], f*Z + g*Y, {f, g}, {Y, Z}]
(*search backward*)NDSolve[{X'[u] == a*X[u] + b Y[u] /. ab, 
Y'[u] == c*Y[u] + d*X[u] + e*Z[u] /. cde, 
Z'[u] == f*Z[u] + g *Y[u] /. fg, X[0.0001] == 69721851.824697, Y[0.0001] == 1710041.1571808, Z[0.0001] == 1249648.07662318}, {X[0], Y[0], Z[0]}, {u, 0, 10.01}, MaxSteps -> 20000]
(*find solution of x[t],y[t],z[t] with NDSolveValue*)
sol = NDSolve[{X'[u] == a*X[u] + b Y[u] /. ab, 
Y'[u] == c*Y[u] + d*X[u] + e*Z[u] /. cde, 
Z'[u] == f*Z[u] + g *Y[u] /. fg, X[0.0001] == 69721851.824697, 
Y[0.0001] == 1710041.1571808, 
Z[0.0001] == 1249648.07662318}, {X[0], Y[0], Z[0]}, {u, 0, 10.01}, MaxSteps -> 20000]
(*plot the results*)
Show[{ListPlot[{xdata, ydata, zdata}, AxesLabel -> {"time", "value"}, PlotLegends -> {"xdata", "ydata", "zdata"}, PlotStyle -> {Red, Green, Pink}], ListLinePlot[{Table[sol[[1]], {u, 0.0, 10, 0.01}], Table[sol[[2]], {u, 0, 10, 0.01}], Table[sol[[3]], {u, 0, 10, 0.01}]}, PlotStyle -> {Orange, Blue, Yellow}, PlotLegends -> {"solution.x", "solution.y", "solution.z"}]}]
share|improve this question

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Browse other questions tagged or ask your own question.