Advanced Search

Question: How to plot a result from Diff Eq

Question:How to plot a result from Diff Eq

Rosen Mitrev 16 Maple
differential-equation plotting
June 04 2008
0

Dear colleagues, i am new in Maple. I am solving a system of four differential equations.I need :

1) to plot the second derivative of the theta[1]. I can plot the theta[1] and its first derivative;

2)to plot V40:

Thanks in advance:)

The code is:

restart;

A[1] := Matrix(4, 4, {(1, 1) = cos(theta[1]), (1, 2) = -sin(theta[1]), (1, 3) = 0, (1, 4) = 0, (2, 1) = sin(theta[1]), (2, 2) = cos(theta[1]), (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}): A[12] := Matrix(4, 4, {(1, 1) = cos(-theta[2]), (1, 2) = 0, (1, 3) = sin(-theta[2]), (1, 4) = 0, (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = -sin(-theta[2]), (3, 2) = 0, (3, 3) = cos(-theta[2]), (3, 4) = L[1], (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}): A[23] := Matrix(4, 4, {(1, 1) = cos(theta[3]), (1, 2) = 0, (1, 3) = sin(theta[3]), (1, 4) = L[2], (2, 1) = 0, (2, 2) = 1, (2, 3) = 0, (2, 4) = 0, (3, 1) = -sin(theta[3]), (3, 2) = 0, (3, 3) = cos(theta[3]), (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}): A[34] := Matrix(4, 4, {(1, 1) = cos(theta[4]), (1, 2) = -sin(theta[4]), (1, 3) = 0, (1, 4) = 0, (2, 1) = sin(theta[4]), (2, 2) = cos(theta[4]), (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = 1, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 1}):

with(LinearAlgebra):

with(plots):

A[2] := simplify(A[1].A[12]): A[4] := simplify(A[1].A[12].A[23].A[34]):

r22 := Vector(4, {(1) = (1/2)*L[2], (2) = 0, (3) = 0, (4) = 1}): r44 := Vector(4, {(1) = L[3], (2) = 0, (3) = 0, (4) = 1}):

   

r20 := simplify(Multiply(A[2], r22)): r40 := simplify(Multiply(A[4], r44)):

r20t := subs(theta[1] = theta[1](t), theta[2] = theta[2](t), r20): r40t := subs(theta[1] = theta[1](t), theta[2] = theta[2](t), theta[3] = theta[3](t), theta[4] = theta[4](t), r40):

V20X := diff(r20t[1], t): V40X := diff(r40t[1], t):

V20Y := diff(r20t[2], t): V40Y := diff(r40t[2], t):

V20Z := diff(r20t[3], t): V40Z := diff(r40t[3], t):

V20 := simplify(sqrt(V20X^2+V20Y^2+V20Z^2)): V40 := simplify(sqrt(V40X^2+V40Y^2+V40Z^2)):

Omega011 := Vector(4, {(1) = 0, (2) = 0, (3) = omega[1], (4) = 0}): Omega122 := Vector(4, {(1) = 0, (2) = omega[2], (3) = 0, (4) = 0}):

Omega010 := Multiply(A[1], Omega011): Omega120 := Multiply(A[2], Omega122):

`ΩA010` := Omega010: `ΩA120` := Omega010+Omega120:

`ΩA011` := MatrixInverse(A[1]).`ΩA010`: `ΩA122` := simplify(MatrixInverse(A[2]).`ΩA120`):

J1 := Matrix(4, 4, {(1, 1) = 0, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = 0, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = J1z, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0}): J2 := Matrix(4, 4, {(1, 1) = J2x, (1, 2) = 0, (1, 3) = 0, (1, 4) = 0, (2, 1) = 0, (2, 2) = J2y, (2, 3) = 0, (2, 4) = 0, (3, 1) = 0, (3, 2) = 0, (3, 3) = J2z, (3, 4) = 0, (4, 1) = 0, (4, 2) = 0, (4, 3) = 0, (4, 4) = 0}):

KR := (1/2)*Transpose(`ΩA011`).J1.`ΩA011`+(1/2)*Transpose(`ΩA122`).J2.`ΩA122`:

KT := (1/2)*m[2]*V20^2+(1/2)*m[4]*V40^2:

KRt := subs(theta[2] = theta[2](t), omega[1] = diff(theta[1](t), t), omega[2] = diff(theta[2](t), t), KR):

K := KRt+KT:

P := m[2]*g*r20t[3]+m[4]*g*r40t[3]:

U := K-P:

First DU

DUdomega01 := subs(diff(theta[1](t), t) = omega[1], U):

Second DU

DUdomega02 := subs(diff(theta[2](t), t) = omega[2], U):

Third DU

DUdomega03 := subs(diff(theta[3](t), t) = omega[3], U):

DUdomega3 := diff(DUdomega03, omega[3]):

DUdomega3t := subs(omega[3] = diff(theta[3](t), t), DUdomega3):

DUdt3 := diff(DUdomega3t, t):

Fourth DU

DUdomega04 := subs(diff(theta[4](t), t) = omega[4], U):

DUdomega4 := diff(DUdomega04, omega[4]):

 

DUdomega1 := diff(DUdomega01, omega[1]):

DUdomega1t := subs(omega[1] = diff(theta[1](t), t), DUdomega1):

DUDt1 := diff(DUdomega1t, t):

DUdtheta11 := subs(theta[1](t) = theta[1], U):

DUdtheta1 := diff(DUdtheta11, theta[1]):

de1 := DUDt1-DUdtheta1+c[1]*theta[1](t)+b[1]*(diff(theta[1](t), t)) = 0:

DUdomega2 := diff(DUdomega02, omega[2]):

DUdomega2t := subs(omega[2] = diff(theta[2](t), t), DUdomega2):

DUdt2 := diff(DUdomega2t, t):

DUdtheta12 := subs(theta[2](t) = theta[2], U):

DUdtheta2 := diff(DUdtheta12, theta[2]):

DUdtheta22 := subs(theta[2] = theta[2](t), DUdtheta2):

de2 := DUdt2-DUdtheta22+b[2]*(diff(theta[2](t), t))+c[2]*theta[2](t) = 0:

DUdtheta13 := subs(theta[3](t) = theta[3], U):

DUdtheta3 := diff(DUdtheta13, theta[3]):

DUdtheta33 := subs(theta[3] = theta[3](t), DUdtheta3):

de3 := DUdt3-DUdtheta33 = 0:

DUdomega4t := subs(omega[4] = diff(theta[4](t), t), DUdomega4):

DUdt4 := diff(DUdomega4t, t):

DUdtheta14 := subs(theta[4](t) = theta[4], U):

DUdtheta4 := diff(DUdtheta14, theta[4]):

DUdtheta44 := subs(theta[4] = theta[4](t), DUdtheta4):

de4 := DUdt4-DUdtheta44 = 0:

ic := theta[1](0) = .1, (D(theta[1]))(0) = 0, theta[2](0) = 0, (D(theta[2]))(0) = 0, theta[3](0) = (1/2)*Pi, (D(theta[3]))(0) = .1, theta[4](0) = 0, (D(theta[4]))(0) = 0:

m[2] := 300: m[4] := 600: L[2] := 3: L[3] := 2: g := 9.81: J1z := 1: J2x := 1: J2y := 1: J2z := 1: b[1] := 5500: b[2] := 5500: c[1] := 500000: c[2] := 500000:

res := dsolve({de1, de2, de3, de4, ic}, {theta[1](t), theta[2](t), theta[3](t), theta[4](t)}, numeric, range = 0 .. 2, output = listprocedure):

plots[odeplot](res, [[t, theta[1](t), colour = blue], [t, theta[2](t), colour = black], [t, theta[4](t), colour = green]], t = 0 .. 2, refine = 1):

V40:

Please Wait...