Group Abstract Group Abstract

Message Boards Message Boards

0
|
8.3K Views
|
20 Replies
|
12 Total Likes
View groups...
Share
Share this post:

Solve a system of differential equations?

Posted 7 years ago

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.

POSTED BY: Torebek Zhumabek
20 Replies
POSTED BY: Hans Dolhaine

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

POSTED BY: Hans Dolhaine

the differential equations system is like Möbius surface.

Ok thanks a lot, I should revise some points

POSTED BY: Torebek Zhumabek
POSTED BY: Hans Dolhaine
Anonymous User
Anonymous User
Posted 7 years ago
POSTED BY: Anonymous User
Anonymous User
Anonymous User
Posted 7 years ago
POSTED BY: Anonymous User
POSTED BY: Daniel Lichtblau
Anonymous User
Anonymous User
Posted 7 years ago
POSTED BY: Anonymous User
Anonymous User
Anonymous User
Posted 7 years ago
POSTED BY: Anonymous User

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].

POSTED BY: Torebek Zhumabek

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]
POSTED BY: Hans Dolhaine

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 :)

POSTED BY: Torebek Zhumabek
POSTED BY: Hans Dolhaine

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] == 0,
   u'[t] == -((1/(h[t]*h[t] + 1/4)^3/2) - 1)*h[t],
   u[0] == 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[1][h ];
up = Derivative[1][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}]
POSTED BY: Hans Dolhaine

Thanks a lot, I am trying to work out it.

Could you show how to plot functions h[t] and u[t] please?

POSTED BY: Torebek Zhumabek

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

POSTED BY: Hans Dolhaine

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

Attachment

Attachments:
POSTED BY: Torebek Zhumabek
POSTED BY: Daniel Lichtblau

Welcome to Wolfram Community! Please make sure you know the rules: https://wolfr.am/READ-1ST

The 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]

enter image description here

POSTED BY: EDITORIAL BOARD

What have you tried? What is unclear about the error message?

POSTED BY: Daniel Lichtblau
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard