NUS-MATH570 Finance (Laplace)

We need to solve this system:

\begin{equation} \begin{cases} \dv{I}{t} = -0.73U + 0.0438 + 0.4 \dv{M}{t} \\ \dv{U}{t} = 0.4I - 0.012 \\ \dv{G}{t} = \dv{M}{t}- I \\ M(t)=0.02\sin (1.15t + \phi) \end{cases} \end{equation}

To be able to work on this, let us create some functions:

variable t, dm = var("t dm") # functions I = function("_I")(t) # _I because i is imaginary U = function("U")(t) G = function("G")(t) # parameter phi = var("phi", latex_name="\phi") # our equations eqns = [ diff(I,t) == -0.73U + 0.0438 + 0.4dm, diff(U,t) == 0.4*I - 0.012, diff(G,t) == dm - I ] eqns desolve(eqns, U, ivar=t, algorithm="fricas").expand() Great, now, we will run the laplace transform upon these equations:

laplace variable s = var("s") # laplaced functions Fi = var("Fi") Fu = var("Fu") Fg = var("Fg") Fm = var("Fm") # constants I0, U0, G0, M0 = var("I0 U0 G0 M0") # substitution dictionary subs = { laplace(I,t,s): Fi, laplace(U,t,s): Fu, laplace(G,t,s): Fg, laplace(M,t,s): Fm, I(0): I0, G(0): G0, U(0): U0, M(0): M0, } # laplace eqns laplace_eqns = [i.laplace(t, s).subs(subs) for i in eqns] laplace_eqns <ipython-input-236-2a2ddfe91635>:20: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. I(Integer(0)): I0, <ipython-input-236-2a2ddfe91635>:21: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. G(Integer(0)): G0, <ipython-input-236-2a2ddfe91635>:22: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. U(Integer(0)): U0, <ipython-input-236-2a2ddfe91635>:23: DeprecationWarning: Substitution using function-call syntax and unnamed arguments is deprecated and will be removed from a future release of Sage; you can use named arguments instead, like EXPR(x=..., y=...) See http://trac.sagemath.org/5930 for details. M(Integer(0)): M0, [Fis - I0 == 0.4Fms - 0.73Fu - 0.4M0 + 0.0438/s, Fus - U0 == 0.4Fi - 0.012/s, Fgs - G0 == Fms - Fi - M0, Fm == (0.02ssin(phi) + 0.023cos(phi))/(s^2 + 1.3225)] And then, let us solve the Laplace solutions:

substitute laplace_solutions = solve(laplace_eqns, Fi, Fu, Fg, Fm, solution_dict=True)[0] laplace_solutions {Fi: 1/100(80000(125I0 - 50M0 + sin(phi))s^4 - 2000(3650U0 - 46cos(phi) - 219)s^3 + 200(66125I0 - 26450M0 + 438)s^2 - 193085(50U0 - 3)s + 115851)/(100000s^5 + 161450s^3 + 38617s), Fu: 1/50(5000000U0s^4 + 4000(500I0 - 200M0 + 4sin(phi) - 15)s^3 + 100(66125U0 + 184cos(phi) + 876)s^2 + 26450(100I0 - 40M0 - 3)s + 115851)/(100000s^5 + 161450s^3 + 38617s), Fg: 1/100(200000(50G0 - 50M0 + sin(phi))s^5 - 10000(1000I0 - 400M0 - 23cos(phi) + 8sin(phi))s^4 + 200(80725G0 - 80725M0 + 36500U0 - 460cos(phi) + 292sin(phi) - 2190)s^3 - 40(330625I0 - 132250M0 - 1679cos(phi) + 2190)s^2 + 193085(20G0 - 20M0 + 50U0 - 3)s - 115851)/(100000s^6 + 161450s^4 + 38617s^2), Fm: 2/5(20ssin(phi) + 23cos(phi))/(400s^2 + 529)} Now we inverse Laplace transform:

I_s(t) = inverse_laplace(laplace_solutions[Fi], s, t) U_s(t) = inverse_laplace(laplace_solutions[Fu], s, t) G_s(t) = inverse_laplace(laplace_solutions[Fg], s, t) M_s(t) = inverse_laplace(laplace_solutions[Fm], s, t) (I_s,U_s,G_s,M_s) (t |--> -1/2061000sqrt(730)(103050U0 + 368cos(phi) - 6183)sin(1/50sqrt(730)t) + 1/1030500(1030500I0 - 412200M0 - 2336sin(phi) - 30915)cos(1/50sqrt(730)t) + 529/51525cos(23/20t)sin(phi) + 529/51525cos(phi)sin(23/20t) + 3/100, t |--> 1/37613250sqrt(730)(1030500I0 - 412200M0 - 2336sin(phi) - 30915)sin(1/50sqrt(730)t) + 1/103050(103050U0 + 368cos(phi) - 6183)cos(1/50sqrt(730)t) - 184/51525cos(phi)cos(23/20t) + 184/51525sin(phi)sin(23/20t) + 3/50, t |--> -1/15045300sqrt(730)(1030500I0 - 412200M0 - 2336sin(phi) - 30915)sin(1/50sqrt(730)t) - 1/41220(103050U0 + 368cos(phi) - 6183)cos(1/50sqrt(730)t) + 1/103050(920cos(phi) + 2061sin(phi))cos(23/20t) + 1/103050(2061cos(phi) - 920sin(phi))sin(23/20t) + G0 - M0 + 5/2U0 - 3/100t - 3/20, t |--> 1/50cos(23/20t)sin(phi) + 1/50cos(phi)sin(23/20t)) Some plots.
I_specific = I_s.subs(I0=0.024, U0=0.039, M0=0, G0=0, phi=2.35) U_specific = U_s.subs(I0=0.024, U0=0.039, M0=0, G0=0, phi=2.35) G_specific = G_s.subs(I0=0.024, U0=0.039, M0=0, G0=0, phi=2.35) M_specific = M_s.subs(I0=0.024, U0=0.039, M0=0, G0=0, phi=2.35) plot(I_specific, t, 0, 10, color="blue") + plot(U_specific, t, 0, 10, color="orange") + plot(G_specific, t, 0, 10, color="green") + plot(M_specific, t, 0, 10, color="red") /Users/houliu/.sage/temp/baboon.jemoka.com/16964/tmp_sei9raar.png