Re: NDSolve and differential equation system
- To: mathgroup at smc.vnet.net
- Subject: [mg54027] Re: NDSolve and differential equation system
- From: "Jens-Peer Kuska" <kuska at informatik.uni-leipzig.de>
- Date: Tue, 8 Feb 2005 05:31:02 -0500 (EST)
- Organization: Uni Leipzig
- References: <cu1vjh$ltj$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
Hi, NDSolve[] has a step size control, any numerical integration assume that the right hand side has smooth derivatives and use this information to estimate the step size. Your UnitStep[] functions introduce nonsmooth derivatives at the right hand side and mix up the step size control. You can a) use the wonderful advanced interface of NDSolve[] to use a constant step size or b) use smooth functions like Erf[] or ArcTan[] instead of UnitStep[] Regards Jens "Christian Moosmann" <some@Adress> schrieb im Newsbeitrag news:cu1vjh$ltj$1 at smc.vnet.net... > 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 >