# Solve a system of differential equations?

Posted 2 months ago
1052 Views
|
20 Replies
|
12 Total Likes
|
 Hello. I need to solve differential equations system dh/dt=f(h) d^2h/dt^2=g(h), h=h(t) Here is the code: ksi = 1; system = {h'[t] == ksi/2*((1/(h[t]*h[t] + 1/4)^3/2) - 1), h''[t] == -((1/(h[t]*h[t] + 1/4)^3/2) - 1)*h[t]}; sol = DSolve[system, {h[t]}, t]; Plot[Evaluate[{h[t]} /. sol], {t, -1, 10}, WorkingPrecision -> 20] And it shows an error. Could anybody fix this problem? Thanks a lot. Answer
20 Replies
Sort By:
Posted 2 months ago
 What have you tried? What is unclear about the error message? Answer
Posted 2 months ago
 Welcome to Wolfram Community! Please make sure you know the rules: https://wolfr.am/READ-1STThe rules explain how to format your code properly. If you do not format code, it may become corrupted and useless to other members. Please EDIT your post and make sure code blocks start on a new paragraph and look framed and colored like this. int = Integrate[1/(x^3 - 1), x]; Map[Framed, int, Infinity]  Answer
Posted 2 months ago
 ksi = 1; system = {h'[t] == ksi/2((1/(h[t]h[t] + 1/4)^3/2) - 1), h''[t] == -((1/(h[t]h[t] + 1/4)^3/2) - 1)h[t]}; sol = DSolve[system, {h, h'}, t]; Plot[[{h[t]} /. sol], {t, -1, 10}, WorkingPrecision -> 20]I am new to Mathematica, and it shows the following error message  Attachments: Answer
Posted 2 months ago
 I will suggest a few things.(1) Clear all Global context symbols (or restart the kernel).(2) Remove the second derivative equation. It makes the system overdetermined.(3) If you want to plot a solution, give an initial value in the system.Beyond that, it is possible the system has no (known) exact solution, so be prepared to use NDSolve instead. Answer
Posted 2 months ago
 Seems your system is incosistent:With hprime = 1/2 (-1 + 1/(2 (1/4 + h[t]^2)^3)); hdoubleprime = h[t] (1 - 1/(2 (1/4 + h[t]^2)^3)); and FullSimplify[Together[D[hprime, t] - hdoubleprime]] I don't get zero Answer
Posted 2 months ago
 But you can do the following. Your system is incosistent, but you could ask for another function u. Redefine your system and solve it numerically (NDSolve) ksi = 1; system = { h'[t] == ksi/2*((1/(h[t]*h[t] + 1/4)^3/2) - 1), h == 0, u'[t] == -((1/(h[t]*h[t] + 1/4)^3/2) - 1)*h[t], u == 0 }; sol = NDSolve[system, {h, u}, {t, 0, 10}] // Flatten; Get the functions h and u h = h /. sol; u = u /. sol; And their derivatives hp = Derivative[h ]; up = Derivative[u]; Calculate values of these at certain values of t hppts = Table[{t, hp[t]}, {t, 0, 2, .05}]; uppts = Table[{t, up[t]}, {t, 0, .6, .05}]; and look whether these values fit to the right-hand sides of your differential equation(s) Plot[ksi/2*((1/(h[t]*h[t] + 1/4)^3/2) - 1), {t, 0, 1}, PlotRange ->All, Epilog -> {Red, PointSize[.015], Point /@ hppts}] and Plot[-((1/(h[t]*h[t] + 1/4)^3/2) - 1)*h[t], {t, 0, .6}, PlotRange -> All, Epilog -> {Red, PointSize[.015], Point /@ uppts}] Answer
Posted 2 months ago
 Thanks a lot, I am trying to work out it.Could you show how to plot functions h[t] and u[t] please? Answer
Posted 2 months ago
 Seems to be a chemical system? Note that in "system" you may choose other values for h[ 0 ] and u[ 0 ]. Concerning the plots try Plot[{h[t], u[t]}, {t, 0, 1.5}, PlotRange -> All] Answer
Posted 2 months ago
 Thanks a lot again. Actually I am doing PhD in celestial mechanics, and this is dynamical equations of the restricted three-body problem. I am From Almaty, Kazakhstan.I need to solve the differential equations system and obtain plots of h[t] and u[t]. Trying to work out :) Answer
Posted 2 months ago
 OK. But what do you expect? Are you sure your system (which is inconsistent as pointed out above) is correct?Do you really mean h'[t] == ksi/2*((1/(h[t]*h[t] + 1/4)^3/2) - 1) or perhaps h'[t] == ksi/2*((1/(h[t]*h[t] + 1/4)^(3/2)) - 1) How did you arrive at your system of equations? What is the underlying physics? Answer
Posted 2 months ago
 As I consider the restricted three-body problem, three bodies form an equilateral triangle. Variable h(t) is height of this equilateral triangle and massless body is on vertex of this triangle. The differential equations system describes the dynamics of the restricted three-body problem.The system is inconsistent and correct. I expect to plot numerical values of functions h[t] and u[t]. Answer
Posted 2 months ago
 A algebraic linear system is consistent if and only if the rightmost column of the augmented matrix is not a pivot column, that is, if and only if an echelon form of the augmented matrix nhas no row of the form: [ 0 ... 0 b ] (with b non-zero) If a linear system is consistent, then the solution set contains either (i) a uniqure solution, when there are no free variables, or (ii) infinitely man solutions, when there is at least one free variable.As to differential equations: those conditions hold true but are not sufficient to say a differential equation pair is solvable or then if solved representable as elementary functions.An example of an inconsistent by the above is: x=5 x=4 Which has no solution possible.(the poster above has used terminology differently - expressing something that is inconsistent could have solutions) Answer
Posted 2 months ago
 i have not reviewed the ODE properties of the 2 equations i'm not swift at reading them in mathematica yeti suggest if there is some discussion over validity that the original poster plots (isoclines) of the two (differential equations) and see visually (though it's possible to be "tricked" by such plots)https://reference.wolfram.com/language/howto/PlotTheResultsOfNDSolve.htmlEquationTrekkeri would say if Hans Dolhaine said there is a problem then to believe him Answer
Posted 2 months ago
 "The system is inconsistent and correct."A more accurate statement would be "the system is inconsistent and (therefore) unsolvable". So this model needs to be rethought. Answer
Posted 2 months ago
 the first problem (for me) is the y''. i know it's possible to reduce an any-order system to a first level system but do not explain or show proof here.substitution may be helpfulif i let ksi=1 u==(1/(y^2 + 1/4))^(3/2) y' == -1+1/2u y'' == (1 - 1/u) y the first diffeq is a linear non-homogenous first order ordinary differential equation linear in u, F(u,y')=0 and F(u,y'')=Q(x) : when solved u will have to be substituted and resolved.the second is a second order equation of the same properties but it is homogenous (and needs u substituted afterward).i haven't tried solving either of the above: they may be too complicated to solve by hand or non-linear after u is found, idkonce either is solved it could be substituted for in the other equation by a grab bag of means. as y (with precaution), finding y' one can find from it y'' and substitute that. doing so might make the system easier to understandthat's as far as i go - i got homework to do :)i can say if you solve a system of linear equations using series method then you'd prefer having initial conditions that method needs, where 0 suggests maclaurin series, and that ultimately NDsolve and numeric approximation should be understood because it can do allot for you that DSolve cannot: if approximation methods understood, if not NDSolve might serve to confuse. Answer
Posted 2 months ago
 A for what Hans suggested I hadn't gotten as far and didn't mention ...definition and solutions of a system of first order equations. 31.1. the pair of equations: dy1/dt = f1(y1,y2,...,yn,t) = f1(t)x+g1(t)y+h1(t) dy2/dt = f2(y1,y2,...,yn,t) = f2(t)x+g2(t)y+h2(t) ... where f1 and f2 are functions of x,y,t fefined on a common set S, is called a system of two first order equations, the sol'n will be two functions on common interval I contained in S satisfying both.means you can solve dx/dt=t/x^2 , dy/dt=y/t^2 does not mean you are encouraged to write (for the sake of a book chapter) dy/dx=y+2y d^2y/dx^2=2y+2Y because that is dy1/dt = f1(y1,t) = f1(t)y1+h1(t) dy1/dt = f2(y1,t) = f2(t)y1+h2(t) while the following may be an identity in x, it may be (is in most cases) inconsistent x' = f1(x) x' = f2(x) also a note that in "beginner's ODE" it is spelled out that fn(t) has no exponent restriction and must be the independent variable tthat being said i see why Hans decided y' and y'' are likely not the same.i question Hans a little as since you said it's not a copied equation but your own (observation): maybe Mathematica's answer is the answer you seek. but i have not bench checked it.you said your studying celestial mechanics, perhaps diff equation course is prerequisite (if it is not then you shouldn't need the answer the question), perhaps you did take the course a while back and need help remembering?it's a fairly complex topic so unless the answer is a single function (often not), describing the answer would be difficult unless you had taken the course. solution, unless otherwise stated, means a complete and exact function but the answer may be implicit or have hidden answers and know these is the beginning of the course. Answer
Posted 2 months ago
 I thought about your problem and - may be I am wrong - would sayIf you have two huge masses M positioned at { - x, 0 } and { x, 0 } and if these masses are effectifely immobile and if you have a tiny mass m at { 0, h } these three form an isoceles triangle and I get as equation of motion for m h''[ t ] == - G M h[ t ] / Sqrt[ x^2 + h[ t ]^2 ]^3 So, where is your 2nd equation or what do you mean by it ?The system described before is essentially sort of a pendulum: the tiny mass oscillates between the two huge masses: Clear[h] sol = NDSolve[{h''[t] == -h[t]/Sqrt[1 + h[t]^2]^3, h == 1, h' == 0}, h, {t, 0, 10}] h = h /. Flatten[sol] Plot[h[t], {t, 0, 10}] Answer
Posted 2 months ago
 the differential equations system is like Möbius surface.Ok thanks a lot, I should revise some points Answer
Posted 2 months ago
 The "pendulum" described above is sort of unphysical. The two huge masses can't stay at rest, they should be attracted by gravitation. So here is another system, three masses (restricted to a plane) but following Newtons law n = 3;(* number of masses *) dim = 2; (* dimension of space - here a plane *) pos = Table[x[i, j][t], {i, 1, n}, {j, 1, dim}]; mass = {10^-4, 1., 1.05}; (* Params of the system *) grav = 1.; (* constant of gravity *) force[i_] := - Sum[If[j == i,0, (grav mass[[j]])/ Sqrt[(pos[[i]] - pos[[j]]).(pos[[i]] - pos[[j]])] (pos[[i]] - pos[[j]])], {j, 1, n}] Initial conditions initpos = Thread /@ Thread[(pos /. t -> 0) == {{0, 5}, {-2, 0}, {2, 0}}]; initspeed = Thread /@ Thread[(D[pos, t] /. t -> 0) == {{0, .2}, {0, -.15}, {0, .15}}]; ODE and its solution sys = Flatten[Join[DGL, initpos, initspeed]]; sol = Flatten[NDSolve[sys, Flatten[pos /. a_[t] -> a], {t, 0, 200}]]; xx = Partition[(Flatten[pos]) /. sol, dim]; and look at it Animate[ Graphics[{PointSize[.05], Point /@ xx /. t -> tt}, PlotRange -> {{-6, 6}, {-10, 10}}], {tt, 0, 50}] It seems that the small mass finally picks up so much speed that it leaves the system Answer
Posted 2 months ago
 The equation for force[ i ] is not correct. It must read force[i_] := -Sum[ If[j == i, 0, (grav mass[[j]])/ Sqrt[(pos[[i]] - pos[[j]]).(pos[[i]] - pos[[j]])] ^3 (pos[[i]] - pos[[j]])], {j, 1, n}] Answer
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.