# DSolve fails to solve a very simple boundary value problem

Posted 10 years ago
9274 Views
|
10 Replies
|
1 Total Likes
|
 Hi there. I was surprised by the following issue with DSolve[]. I solve a simple boundary value problem:A v''(t)+B v(t) = F(t)where A, B -- 3x3 constant matrices, v(t) -- 3-component vector-function to be search, F(t) -- given vector-function. Each component of v(t) satisfies boundary conditions at the point 0 and 1. The problem is straightforward and it seems no issue of using DSolve should appear for any admissible F(t).In my case is calculated on previous steps. The right part comprises polynomials and hyperbolic sine/cosine.A={{-0.044697669173699448530657096189269774609066666666666666666667437866249.52287874528033,-0.012243540763682102904177110474279025634133333333333333333333563124849.52287874528033,0},{-0.012243540763682102904177110474279025634133333333333333333333563124849.52287874528033,-0.053659814734196655264138211329287943019266666666666666666667592496249.52287874528033,-0.009622111232315166943183972182021976285866666666666666666666847258249.52287874528033},{0,-0.009622111232315166943183972182021976285866666666666666666666847258249.52287874528033,-0.033141781673581725468679476493387322164266666666666666666667238484549.52287874528033}};B={{0.134093007521098345591971288567809323827200000000000000000002381302349.52287874528033,-0.073461244582092617425062662845674153804800000000000000000001304567849.52287874528033,0},{-0.073461244582092617425062662845674153804800000000000000000001304567849.52287874528033,0.160979444202589965792414633987863829057800000000000000000002858767349.52287874528033,-0.057732667393891001659103833092131857715200000000000000000001025250649.52287874528033},{0,-0.057732667393891001659103833092131857715200000000000000000001025250649.52287874528033,0.099425345020745176406038429480161966492800000000000000000001765653549.52287874528033}};Cv={0.06063176293900572816690862572213517002239999999959703682741.3422565562137*t,0.21278032964063766549749999999999999999999999999992032228440.30310235657392-0.03321638353144512597695900434870819654034999999987569840340.30279030542798*t+0.0160691980909371055151909786571919738642539.93202436395856*Cosh[1.73205080756887729352744634150587236694280525381038062805647.307206068870755*t]-0.03776972085899651059900806241041891399654810098499796775840.303173934606086*Sinh[1.73205080756887729352744634150587236694280525381038062805647.307206068870755*t],0.111481694136719805460000000000000000000000000000000000000001979757349.69897000433602};(here Cv used instead of F(t) ). For solving the problem the following code is used:Q = DSolve[{A[[1,1]] v1''[t]+A[[1,2]] v2''[t]+B[[1,1]] v1[t]+B[[1,2]] v2[t]==Cv[[1]],  A[[2,1]] v1''[t]+A[[2,2]] v2''[t]+A[[2,3]] v3''[t]+B[[2,1]] v1[t]+B[[2,2]] v2[t]+B[[2,3]] v3[t]== Cv[[2]],A[[3,2]] v2''[t]+A[[3,3]] v3''[t]+B[[3,2]] v2[t]+B[[3,3]] v3[t]== Cv[[3]], v1[0]  ==  0, v1[1] == 1, v2[0]  ==  0, v2[1] == 1, v3[0]  ==  0, v3[1] == 1}, {v1[t],v2[t],v3[t]}, t]But calculation invokes the errorNo more memory available.Mathematica kernel has shut down.Try quitting other applications and then retry.or the error of kernelI tried to transform the right part into exponential formCv0=TrigToExp[Cv]Now.But DSolve cannot solve the problem at all:Q = DSolve[{A[[1,1]] v1''[t]+A[[1,2]] v2''[t]+B[[1,1]] v1[t]+B[[1,2]] v2[t]== Cv[[1]], A[[2,1]] v1''[t]+A[[2,2]] v2''[t]+A[[2,3]] v3''[t]+B[[2,1]] v1[t]+B[[2,2]] v2[t]+B[[2,3]] v3[t]== Cv0[[2]],A[[3,2]] v2''[t]+A[[3,3]] v3''[t]+B[[3,2]] v2[t]+B[[3,3]] v3[t]== Cv0[[3]], v1[0]  ==  0, v1[1] == 1, v2[0]  ==  0, v2[1] == 1, v3[0]  ==  0, v3[1] == 1}, {v1[t],v2[t],v3[t]}, t]Result:Indeed, I don't understand, what goes wrong?
10 Replies
Sort By:
Posted 10 years ago
 This produced an answer after several minutes and using no more than 400MB of memory.  I tried solving the BVP directly but it took more than an hour longer; produced warning/error messages about exponents being too large; and it has an answer more than 20 times as large and takes much longer to evaluate.  So I opt to present this way. sol0 = DSolve[    SetPrecision[{A.{x''[t], y''[t], z''[t]} + B.{x[t], y[t], z[t]} ==        Cv}, Infinity],    {x, y, z}, t];   Block[{Integrate = NIntegrate, eqns},  eqns = {x[0] == 0, x[1] == 1, y[0] == 0, y[1] == 1, z[0] == 0, z[1] == 1} /. sol0;  NSolve[eqns, Cases[eqns, _C, Infinity] // Union] ](*  {{C[1] -> -0.0249078, C[2] -> 1.64911, C[3] -> 1.3804,     C[4] -> -1.22953, C[5] -> 1.48102, C[6] -> -1.06965}}*)Taking a hint from Sean Clarke's observation, I changed the coefficients to exact numbers.  Using  SetPrecision[#, Infinity]  is not always a good idea.  All the real numbers get converted to exact fractions with very large denominators.  I first tried i out the idea by rounding the coefficients to the nearest 1/1000.  That worked, but the answer is in terms of Integrate and RootSum objects.  It made me think this problem is not as simple as some have supposed.  It is also why I used  Block  to switch  Integrate  to  NIntegrate.After getting the general solution, we can solve for the constants that satisfy the boundary conditions.If you need 40-50 digits of precision as the OP's input suggests, then the above will have to be modified.Finally, I would think Bill Simpson's approach is best, unless one needs to analyze the algebraic expressions in a formula that occupies 300-400K of memory.
Posted 10 years ago
 Konstantin, right, sorry, I've meant first order.Sure, you can  easily transform you system into first order one.I.M.
Posted 10 years ago
 Certainly, Ivan...
Posted 10 years ago
 A = <<>>; B = <<>>; Cv =<<>>; sol = {v1[t], v2[t], v3[t]} /.     NDSolve[{A[[1, 1]] v1''[t] + A[[1, 2]] v2''[t] + B[[1, 1]] v1[t] + B[[1, 2]] v2[t] == Cv[[1]],      A[[2, 1]] v1''[t] + A[[2, 2]] v2''[t] + A[[2, 3]] v3''[t] + B[[2, 1]] v1[t] + B[[2, 2]] v2[t] + B[[2, 3]] v3[t] == Cv[[2]],      A[[3, 2]] v2''[t] + A[[3, 3]] v3''[t] + B[[3, 2]] v2[t] + B[[3, 3]] v3[t] == Cv[[3]], v1[0] == 0, v1[1] == 1,      v2[0] == 0, v2[1] == 1, v3[0] == 0, v3[1] == 1}, {v1[t], v2[t], v3[t]}, {t, 0, 1}][[1]]; Plot[sol, {t, 0, 1}]and in a fraction of a second the three functions are plotted.Now how accurate are the three plots? At least the end points appear to be correct.
Posted 10 years ago
 That is indeed a simple problem. I strongy suspect that DSolve[] has nothing to do with it, it's just not enough physical memory on your system to deal with this task. Try to monitor your system resources while varying he precision of operations.Anyway, the linear system of ODEs, can be solved without any problems using DSolve[] with symbolic coefficients, then numbers can be plugged  into obtained solution.Moreover, there is no need to use DSolve[] at all. The solutuon for linear system of the formz(t)'' + A z(t) = F(t) (original)z(t)' + A z(t) = F(t)  (corrected)is justz(t) = Exp(-A t) z(0) + Exp(-A t) Int[( f(x) Exp(A x), {x,0,t}](not sure about the signs in the integral), where A is a matrix.The above is also true for non-autonomous systems with A(t), where it is only necessary to extend the system.I.M.
Posted 10 years ago
 Ivan, you are wrong in this case, because my ODE is of 2nd order:But this approach can be used for 1st order ODE. My ODE can easily be transformed into such a ODE and solved. Nevertheless, I would like to use DSolve, because it requires much less job. If DSolve or NDSolve will demonstrate full frustration, I'll turn to other methods.
Posted 10 years ago
 I see, Sean. It is very pity. Such powerful CAS cannot solve a simple, in principle, problem. True misfortune.
Posted 10 years ago
 It'll still suffer from the same fundamental problems.  I honestly don't try to guess too much about the internal workings of DSolve. DSolve will very often flatly refuse to start working with floating point numbers. Maybe it figures that exponentials and floating point numbers are a potentially dangerous coupling and so refuses to try anything.
Posted 10 years ago
 Thanks, Sean. I suspected the points you listed. But can you suggest, why DSolve doesn't work in my last example, where the right part was converted to exponential form? Maybe, the expression has become very complicated for DSolve?
Posted 10 years ago
 DSolve very often doesn't like to use Floating point numbers. This is because it is a symbolic function and may not perform numerically stable operations on the floating point numbers in the course of its symbolic operation. Be careful when mixing symbolic and numeric computations. Generally you might consider: 1. Turning the floating point numbers into rational numbers with Rationalize2. Numerically Solving the equation with NDSolve3. Solve using Parameters and then substitute the values of those parameters into the result from DSolve.If you do (3), try running the following:DSolve[{B[1,1] v1[t]+B[1,2] v2[t]+A[1,1] (v1^\[Prime]\[Prime])[t]+A[1,2] (v2^\[Prime]\[Prime])[t]==1,2 v1[t]+2 v2[t]+2 v3[t]+2 (v1^\[Prime]\[Prime])[t]+2 (v2^\[Prime]\[Prime])[t]+2 (v3^\[Prime]\[Prime])[t]==2,3 v2[t]+3 v3[t]+3 (v2^\[Prime]\[Prime])[t]+3 (v3^\[Prime]\[Prime])[t]==3,v1[0]==0,v1[1]==1,v2[0]==0,v2[1]==1,v3[0]==0,v3[1]==1},{v1[t],v2[t],v3[t]},t]This solves quickly and you can decide how to substitute the actual values back in. This can be tricky sometimes because of numerical precision issues.