Re: NDSolve and differential equation system

*To*: mathgroup at smc.vnet.net*Subject*: [mg54039] Re: NDSolve and differential equation system*From*: Maxim <ab_def at prontomail.com>*Date*: Tue, 8 Feb 2005 05:31:28 -0500 (EST)*Organization*: MTU-Intel ISP*References*: <cu1vjh$ltj$1@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

On Sat, 5 Feb 2005 08:20:01 +0000 (UTC), Christian Moosmann <some@Adress> wrote: > Hi MathGroup > > I have a question concerning NDSolve and a system of differential > equations. > > consider an equation system A x'[t] + (K1 + v[t] * K2).x[t]=B, where > A,K1,K2 are Matrices, B and x are vectors and v[t] is a defined skalar > function. I want to compute this system with NDSolve, however, it takes > a very long time. > > To try it you should use Theodor Gray's ShowStatus function in the > frontend. > > ShowStatus[status_] := > LinkWrite[$ParentLink, > SetNotebookStatusLine[FrontEnd`EvaluationNotebook[], > ToString[ status ] ] ]; > > Dim = 40; > > I use random Matrices here. They show the same basic behaviour as the > matrices I usually use (those are also dense), so this problem should > not be related directly to the matrices. > > K1 = Table[Random[], {i, Dim}, {j, Dim}]; > K2 = Table[Random[], {i, Dim}, {j, Dim}]; > Ainv = Inverse[Table[Random[], {i, Dim}, {j, Dim}]]; > B = Table[Random[], {i, Dim}]; > > tStart = 0.; > tEnd = 1; > > At first we solve without the function v[t]: > > Timing[NDSolve[{D[ x[t], t] + Ainv.( K1 + K2).x[t] == B , > x[tStart] == Table[0. , {Dim}]} , x , { t, tStart, tEnd }, > SolveDelayed -> True, StepMonitor :> ShowStatus[t] ] ] > > {0.080987 Second, {{x -> InterpolatingFunction[{{0., 1.}}, <>]}}} > > is quite fast, now we would like to include the function: > > v[t_] := UnitStep[t - 0.2]*4*(t - 0.2) - UnitStep[t - 0.45]*4*(t - 0.45); > > Timing[NDSolve[{D[ x[t], t] + Ainv.( K1 + v[t]*K2).x[t] == B , > x[tStart] == Table[0. , {Dim}]} , x , { t, tStart, tEnd }, > SolveDelayed -> True, StepMonitor :> ShowStatus[t] ] ] > > {137.341 Second, {{x -> InterpolatingFunction[{{0., 1.}}, <>]}}} > > Well, this is a factor of > 1000, and what I would like to compute is at > least Dimension 100. If you use the stepmonitor you also may notice, > that the time counts up, then stops for some time, counts up again... > > > Does anyone have any suggestions how to deal with that? > > Thanks in advance > > Christian > I managed to do it numerically only up to Dim=80 by rewriting the vector ODE as a system of scalar equations and using Method -> IDA: In[1]:= SeedRandom[1] Dim = 80; K1 = Table[Random[], {Dim}, {Dim}]; K2 = Table[Random[], {Dim}, {Dim}]; Ainv = Inverse[Table[Random[], {Dim}, {Dim}]]; B = Table[Random[], {Dim}]; v[t_] := UnitStep[t - 2/10]*4*(t - 2/10) - UnitStep[t - 45/100]*4*(t - 45/100); Module[ {Lf = Array[x, {Dim}], x}, x = Evaluate[Through[Lf[#]]]&; NDSolve[Thread /@ {x'[t] + Ainv.(K1 + v[t]*K2).x[t] == B, x[0] == Table[0, {Dim}]}, Lf, {t, 0, 1}, Method -> IDA, SolveDelayed -> True] ]; // Timing Out[8]= {175.625 Second, Null} In[9]:= MaxMemoryUsed[] Out[9]= 695334616 This might be a case where the piecewise handling routines of Mathematica 5.1 do more harm than good, because if we choose a different v[t], say v[t_] = (Erf[10(t - 13/40)] + 1)/2, which is very close to the original v[t], then NDSolve can handle this equation with Dim=80 much easier -- both time and memory requirements will be significantly lower. However, if Ainv.K1 and Ainv.K2 commute or, equivalently, if K1.Ainv.K2 == K2.Ainv.K1 (Ainv is nonsingular), then the solution can be written out explicitly by analogy with the solution of the scalar equation x'[t] + a[t]*x[t] == b. We still won't be able to find the matrix exponentials symbolically, but we can use them to approximate the solution. Let's take K2 = -2*K1, which satisfies the above condition: In[1]:= (*starting with a fresh kernel*) SeedRandom[1] Dim = 80; K1 = Table[Random[], {Dim}, {Dim}]; K2 = -2*K1; Ainv = Inverse[Table[Random[], {Dim}, {Dim}]]; B = Table[Random[], {Dim}]; v[t_] := UnitStep[t - 2/10]*4*(t - 2/10) - UnitStep[t - 45/100]*4*(t - 45/100); Lsol = Module[ {Lf = Array[x, {Dim}], x, V, f, df, g, gv, ig, Fv, dFv}, x = Evaluate[Through[Lf[#]]]&; V[t_] = Integrate[v[tau], {tau, 0, t}, Assumptions -> t > 0]; f[t_] := f[t] = MatrixExp[-Ainv.(K1*t + K2*V[t])]; df[t_] := df[t] = -Ainv.(K1 + K2*v[t]).f[t]; g[t_] := g[t] = MatrixExp[Ainv.(K1*t + K2*V[t])].B; gv[t_?NumericQ, i_] := g[t][[i]]; ig[t_] := ig[t] = NIntegrate[ Evaluate@ Table[gv[tau, i], {i, Dim}], {tau, 0, Min[t, 1/5], Min[t, 9/20], t}]; Fv[t_?NumericQ, i_] := (f[t].ig[t])[[i]]; dFv[t_?NumericQ, i_] := (df[t].ig[t] + f[t].g[t])[[i]]; Table[Interpolation[ Table[{t, {Fv[t, i], dFv[t, i]}}, {t, 0., 1., .01}]], {i, Dim}] ]; // Timing Out[8]= {25.266 Second, Null} In[9]:= MaxMemoryUsed[] Out[9]= 19204344 Some jumping through hoops is required because we cannot declare the arguments of NIntegrate as vectors. An interesting observation is that Interpolation seems to work much better than FunctionInterpolation in this case. This method scales much better, and the agreement between the two solutions is reasonably good: Lsol2 = Module[ {Lf = Array[x, {Dim}], x}, x = Evaluate[Through[Lf[#]]]&; NDSolve[Thread /@ {x'[t] + Ainv.(K1 + v[t]*K2).x[t] == B, x[0] == Table[0, {Dim}]}, Lf, {t, 0, 1}, Method -> IDA, SolveDelayed -> True] ]; Plot[Module[ {Lf = Array[x, {Dim}], x}, x = Evaluate[Through[Lf[#]]]&; Norm[Through[Lsol[t]] - (x[t] /. Lsol2[[1]])] ], {t, 0, 1}, PlotRange -> All] As to why NDSolve issues the NDSolve::ndfdmc warning when the equation is given in the vector form with SolveDelayed set to False, this is most likely due to the incorrect preprocessing of the input: thus, the next two examples have well-defined solutions, but NDSolve gives a warning in the first case and just crashes the kernel in the second case: In[10]:= NDSolve[ {x''[t].x''[t] == 1, x''[t].{-1, 1} == 0, x[0] == {1, 1}, x'[0] == {0, 0}}, x, {t, 0, 1}] NDSolve::overdet: There are fewer dependent variables, {x[t]}, than equations, so the system is underdetermined. Out[10]= NDSolve[{Derivative[2][x][t] . Derivative[2][x][t] == 1, Derivative[2][x][t] . {-1, 1} == 0, x[0] == {1, 1}, Derivative[1][x][0] == {0, 0}}, x, {t, 0, 1}] NDSolve[ {x''[t].x''[t] + x[t].x[t] == 0, x[0] == {1}, x'[0] == {0}}, x, {t, 0, 1}, SolveDelayed -> True] (*crashes the kernel*) It is interesting to compare this with the behaviour of the previous versions of Mathematica: in version 5.0 DSolve could return all kinds of weird output when the equations were given as lists. Now DSolve just generates the warning message DSolve::nolist, but NDSolve has to work with vector equations, and the processing of vector input is still far from perfect: In[11]:= NDSolve[ {{x'[t], y'[t]} == {1, 1}, {x[0], y[0]} == {0, 0}}, {x, y}, {t, 0, 1}] NDSolve::underdet: There are more dependent variables, {x[t], y[t]}, than equations, so the system is underdetermined. Out[11]= NDSolve[{{Derivative[1][x][t], Derivative[1][y][t]} == {1, 1}, {x[0], y[0]} == {0, 0}}, {x, y}, {t, 0, 1}] In the last case the problem can be resolved simply by applying Thread. Maxim Rytin m.r at inbox.ru