- #1
Barnak
- 63
- 0
I need help with the Maple maths software, which I don't know very well. I need to numerically solve a system of three second order differential equations, specifically with Maple. I have to write a code which numerically solves the Newton's equation, using some initial conditions, and then show the numerical trajectory as a 3D plot. How am I supposed to do that?
More specifically, consider the following three coupled differential equations:
[itex]\frac{d^2 x(t)}{dt^2} = F_x(x(t), y(t), z(t))[/itex]
[itex]\frac{d^2 y(t)}{dt^2} = F_y(x(t), y(t), z(t))[/itex]
[itex]\frac{d^2 z(t)}{dt^2} = F_z(x(t), y(t), z(t))[/itex]
where [itex]F_x(x, y, z)[/itex], [itex]F_y(x, y, z)[/itex] and [itex]F_z(x, y, z)[/itex] are three specific functions (which should be specified in the Maple code). The initial conditions are also given:
[itex]x(0) = C_1[/itex],
[itex]y(0) = C_2[/itex],
[itex]z(0) = C_3[/itex],
[itex]v_x(0) = C_4[/itex],
[itex]v_y(0) = C_5[/itex],
[itex]v_z(0) = C_6[/itex].
Here, the six constants [itex]C_i[/itex] are arbitrary, but they should be specified as some real numbers in the Maple code.
How to program Maple so it can solves the previous equations and show the trajecory in some 3D graph ? I need a complete and explicit code to do this, because I'm really a novice with Maple.
In the particular case of Mathematica (a rival of Maple), I know how to do it. The following code gives the complete solution to my problem, except that I need the Maple code which can do the same. How to translate this code to Maple language ? Please, be specific. The answer will certainly help a lot of other people!
To trace the trajectory in 3D, the Mathematica code is bellow :
More specifically, consider the following three coupled differential equations:
[itex]\frac{d^2 x(t)}{dt^2} = F_x(x(t), y(t), z(t))[/itex]
[itex]\frac{d^2 y(t)}{dt^2} = F_y(x(t), y(t), z(t))[/itex]
[itex]\frac{d^2 z(t)}{dt^2} = F_z(x(t), y(t), z(t))[/itex]
where [itex]F_x(x, y, z)[/itex], [itex]F_y(x, y, z)[/itex] and [itex]F_z(x, y, z)[/itex] are three specific functions (which should be specified in the Maple code). The initial conditions are also given:
[itex]x(0) = C_1[/itex],
[itex]y(0) = C_2[/itex],
[itex]z(0) = C_3[/itex],
[itex]v_x(0) = C_4[/itex],
[itex]v_y(0) = C_5[/itex],
[itex]v_z(0) = C_6[/itex].
Here, the six constants [itex]C_i[/itex] are arbitrary, but they should be specified as some real numbers in the Maple code.
How to program Maple so it can solves the previous equations and show the trajecory in some 3D graph ? I need a complete and explicit code to do this, because I'm really a novice with Maple.
![Cry :cry: :cry:](/styles/physicsforums/xenforo/smilies/cry.png)
In the particular case of Mathematica (a rival of Maple), I know how to do it. The following code gives the complete solution to my problem, except that I need the Maple code which can do the same. How to translate this code to Maple language ? Please, be specific. The answer will certainly help a lot of other people!
Code:
Pos[ x_, y_, z_] := {x, y, z}
r[x_, y_, z_] := Sqrt[x^2 + y^2 + z^2]
Velocity[t_] := {x'[t], y'[t], z'[t]}MagnMoment := {0, 0, 1}
MagneticField[ x_, y_, z_] := 3 (MagnMoment.Pos[x, y, z]) Pos[x, y, z]/r[x, y, z]^5 - MagnMoment/r[x, y, z]^5
MagneticCoupling := 0
GravitationCoupling := 1
FrictionCoupling := 0.005
GravitationalForce[t_] := - GravitationCoupling Pos[x[t], y[t], z[t]]/r[x[t], y[t], z[t]]^3
MagneticForce[t_] := MagneticCoupling Cross[Velocity[t], MagneticField[x[t], y[t], z[t]]]
Friction[t_] := - FrictionCoupling Velocity[t]
Force[t_] := GravitationalForce[t] + MagneticForce[t] + Friction[t]
Fx[t_] := {1, 0, 0}.Force[t]
Fy[t_] := {0, 1, 0}.Force[t]
Fz[t_] := {0, 0, 1}.Force[t]
Time1 := 0
Time2 := 200
X0 := 3
Y0 := 0.1
Z0 := 5
Vx0 := -0.5
Vy0 := 0
Vz0 := 0
Trajectory = NDSolve[{
x''[t] == Fx[t],
y''[t] == Fy[t],
z''[t] == Fz[t],
x[0] == X0,
y[0] == Y0,
z[0] == Z0,
x'[0] == Vx0,
y'[0] == Vy0,
z'[0] == Vz0
}, {x,y,z}, {t,Time1,Time2}, Method -> ExplicitRungeKutta, MaxSteps->10000]
To trace the trajectory in 3D, the Mathematica code is bellow :
Code:
PTmin := (x /. Trajectory)[[1]][[1]][[1]][[1]]
PTmax := (x /. Trajectory)[[1]][[1]][[1]][[2]]
ParametricPlot3D[Evaluate[{x[t], y[t], z[t]} /. Trajectory], {t, PTmin, PTmax}, PlotPoints -> 10000]
Last edited: