Numerically solve an IVP

edited August 17 in Assignments

(10 pts)

Produce and IVP of the form
$$
y' = r\,y\,\sin(y), : y(0)=1
$$
by choosing the decimals of $r$ to be the positions in the alphabet of your login name.
For example, my login name is mark and the positions of those letters are 13, 1, 18, and 11. Thus, for me, $r=0.1311811$ and my ODE is
$$
y' = 0.1311811 \, y \, \sin(y).
$$
Find the value of $t$ for which $y(t)=3$.

Comments

  • edited August 19

    So, here's how I might approach it. First, my IVP is

    $$
    y' = 0.1311811 \, y \, \sin(y), \: y(0) = 1.
    $$

    OK, so first I'm going to load the libraries:

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import odeint
    from scipy.interpolate import interp1d
    from scipy.optimize import brentq
    

    Now, let's set up and solve the IVP:

    # Redefine f as a function of both y and t
    def f(y,t):
        return 0.1311811*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Finally, I'll interpolate the solution and solve it:

    y_interp = interp1d(t,y)
    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    
    #Output:
    11.997330524726108
    
  • edited August 20

    ``So my equation (if I use Chris as my name) is: $$ y' = 0.3818919 \, y \, \sin(y), \: y(0) = 1 $$

    My code libraries are:

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import odeint
    from scipy.interpolate import interp1d
    from scipy.optimize import brentq
    

    Then I set up and solved the IVP:

    # Redefine f as a function of both y and t
    def f(y,t):
      return 0.3818919*y*np.sin(y)
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,10,100)`
    # Solve it!
    y = odeint(f, 1, t).flatten()`
    # Plot it!!
    plt.plot(t,y);
    

    Finally, I interpolated the solution:

    y_interp = interp1d(t,y)
    y_interp
    t0 = brentq(lambda y: y_interp(y)-3, 0, 10)
    t0
    

    And got:

    #Output:
    
    4.12140025882665
    
    mark
  • edited August 19

    My IVP was: $$y' = (0.112935)\,y\,\sin(y), : y(0)=1$$

    Loading the libraries:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    
    # A numerical rootfinder
    from scipy.optimize import brentq
    

    Now for solving the IVP:

    # Redefine f as a function of both y and t
    def f(y,t):
    return (0.112935)*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Interpolating the solution:

    y_interp = interp1d(t,y)
    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    #Output: 13.933757759135883
    

    checking graphically:

    plt.plot(t,y)
    plt.plot([t0,t0],[0,6],'k--')
    plt.plot([0,50],[3,3],'k--')
    plt.plot(t0,3,'ro');
    

    mark
  • My equation is:

    $$ y' = 0.1020815135ysin(y) $$

    Here's the code I did:

     # A graphics library
     import matplotlib.pyplot as plt
    
     # Basic numerical tools
     import numpy as np
    
     # Scipy's numerical ODE solver
     from scipy.integrate import odeint
    
     # Convert data generated by odeint to a computable function
     from scipy.interpolate import interp1d
    
     # A numerical rootfinder
     from scipy.optimize import brentq
    
     # Redefine f as a function of both y and t
     def f(y,t):
          return 0.1020815135*y*np.sin(y)
    
      # Some t values that are appropriate to plot the solution
      t = np.linspace(0,20,200)
    
      # Solve it!
      y = odeint(f, 1, t).flatten()
    
      # Plot it!!
      plt.plot(t,y);
    

      odeint(f,1,t)[:8].flatten()
    

    array([1. , 1.00869456, 1.01751339, 1.02645847, 1.03553176,
    1.04473522, 1.0540708 , 1.06354048])

      y_interp = interp1d(t,y)
      y_interp
    
      y_interp(20)
    

    array(3.10782519)

      t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
      t0
    

    15.41536137275634

      y_interp(t0)
    

    array(3.)

      plt.plot(t,y)
      plt.plot([t0,t0],[0,6],'k--')
      plt.plot([0,50],[3,3],'k--')
      plt.plot(t0,3,'ro');
    

    mark
  • edited August 26

    My IVP was

    $$dy/dt=0.42512114y\sin(y), y(0)=1$$

    Using basic, large step-size Euler's method, I got that at $$y(t)=3, t\approx3.62$$

    To get a more precise estimate, I used Mark's code:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    def f(y,t):
        return 0.42512114*y*np.sin(y)
    t = np.linspace(0,8,100)
    y = odeint(f, 1, t).flatten()
    y
    
    plt.plot(t,y)
    

    from scipy.interpolate import interp1d
    from scipy.optimize import brentq
    
    y_func = interp1d(t,y)
    t0 = brentq(lambda t: y_func(t)-3,0,8)
    t0
    

    To find that at $$y(t)=3, t\approx3.702112268685581$$

  • edited August 20

    My IVP using my name Austin is:

    $$y' = 0.1311811 \, y \, \sin(y), : y(0) = 1.$$

    Next I loaded the libraries:

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import odeint
    from scipy.interpolate import interp1d
    from scipy.optimize import brentq
    

    Then I solved the IVP:

    # Redefine f as a function of both y and t
    def f(y,t):
        return 0.1211920914*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Finally I interpolated the solution and solved it:

    y_interp = interp1d(t,y)
    y_interp
    
    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    #Output:
    12.984599840002112
    
    mark
  • edited August 19

    My IVP is:
    $y^{\prime}=.114418523 y \sin (y), y(0)=1$

    And my code and resulting answers are as follows:

    Import necessary libraries:

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import odeint
    from scipy.interpolate import interp1d
    from scipy.optimize import brentq

    def f(y,t):
    return .114418523*y*np.sin(y)
    t = np.linspace(0,20,100)
    y = odeint(f, 1, t).flatten()
    y

    This gives me the array:
    array([1. , 1.01976475, 1.04017289, 1.06124679, 1.08300904,
    1.10548228, 1.12868906, 1.15265169, 1.17739199, 1.20293107,
    1.2292891 , 1.25648492, 1.28453576, 1.31345679, 1.34326076,
    1.37395749, 1.40555336, 1.43805078, 1.47144761, 1.5057366 ,
    1.54090478, 1.57693284, 1.61379466, 1.65145672, 1.68987771,
    1.72900818, 1.7687903 , 1.80915783, 1.85003621, 1.8913428 ,
    1.9329875 , 1.97487337, 2.01689757, 2.05895252, 2.10092721,
    2.14270855, 2.18418305, 2.22523826, 2.26576451, 2.30565631,
    2.34481388, 2.38314434, 2.42056288, 2.45699355, 2.49236994,
    2.52663557, 2.55974407, 2.59165907, 2.62235406, 2.65181188,
    2.68002421, 2.70699093, 2.73271937, 2.75722353, 2.78052327,
    2.80264357, 2.82361368, 2.84346643, 2.86223752, 2.87996487,
    2.89668801, 2.91244763, 2.92728501, 2.94124171, 2.95435914,
    2.96667832, 2.97823957, 2.98908237, 2.99924512, 3.00876505,
    3.01767815, 3.026019 , 3.03382083, 3.04111539, 3.04793301,
    3.0543026 , 3.06025159, 3.06580604, 3.07099064, 3.07582871,
    3.08034231, 3.08455224, 3.08847808, 3.0921383 , 3.09555024,
    3.0987302 , 3.10169349, 3.10445446, 3.1070266 , 3.1094225 ,
    3.11165399, 3.11373212, 3.11566724, 3.11746902, 3.11914651,
    3.12070814, 3.12216183, 3.12351492, 3.12477431, 3.12594641])

    plt.plot(t,y)

    y_func = interp1d(t,y)
    t1 = brentq(lambda t: y_func(t)-1,0,20)
    t1

    0.0

    plt.plot(t,y)
    plt.plot([t1,t1],[0,1.2],'k--')
    plt.plot(t1,1,'ko')

  • edited August 19

    IVP:

    $$
    y\prime = 0.11475121ysin(y)
    $$

    Load the Libraries:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    
    # A numerical rootfinder
    from scipy.optimize import brentq
    

    Set up and solve the IVP:

    # Redefine f as a function of both y and t
    def f(y,t):
    return (0.11475121)*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Interpolate the solution and solve:

    y_interp = interp1d(t,y)
    y_interp
    
    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    
    #Output:
    8.039049419781893
    
    mark
  • edited August 19

    My equation is (using mcespede):

    $$ y'=0.13351916545 y \sin (y) $$

    My library codes are:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    
    # A numerical rootfinder
    from scipy.optimize import brentq
    

    Then I set up and and solved the IVP:

    # Redefine f as a function of both y and t
    def f(y,t):
        return 0.13351916545*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Interpolate the solution and solve it:

    y_interp = interp1d(t,y)
    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    

    # Output:
    11.787222349796691

    mark
  • edited August 19

    My equation is: $y' = (0.11255) * y * sin(y) $

    Next I load the libraries:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    
    # A numerical rootfinder
    from scipy.optimize import brentq
    

    Next setting up and solving the IVP

    # Redefine f as a function of both y and t
    def f(y,t):
    return (0.11255)*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Finally, Interpolating the solution to solve

    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    

    The output being
    13.982315772416888

    mark
  • My equation is: $$y' = 0.51893y\sin(y)$$

    First the following libraries are loaded:

      # A graphics library
      import matplotlib.pyplot as plt
    
      # Basic numerical tools
      import numpy as np
    
      # Scipy's numerical ODE solver
      from scipy.integrate import odeint
    
      # Convert data generated by odeint to a computable function
      from scipy.interpolate import interp1d
    
      # A numerical rootfinder
      from scipy.optimize import brentq
    

    Next the solution is graphed with the following code:

      # Redefine f as a function of both y and t
      def f(y,t):
          return 0.51893*y*np.sin(y)
    
      # Some t values that are appropriate to plot the solution
      t = np.linspace(0,5,1000)
    
      # Solve it!
      y = odeint(f, 1, t).flatten()
    
      # Plot it!!
      plt.plot(t,y);
    

    Finally interpolating the solution using:

      y_interp = interp1d(t,y)
      t0 = brentq(lambda y: y_interp(y)-3, 0, 5)
      t0
    

    Output: 3.0323732170171915

    Displaying this graphically:

    plt.plot(t,y)
    plt.plot([t0,t0],[1,3],'k--')
    plt.plot([0,5],[3,3],'k--')
    plt.plot(t0,3,'ro');
    

    mark
  • edited August 20

    Using my name generates my equation as:
    $$y' = 0.1920516 \, y \, \sin(y), : y(0) = 1$$

    Loading my libraries:

        import matplotlib.pyplot as plt
        import numpy as np
        from scipy.integrate import odeint
        from scipy.interpolate import interp1d
        from scipy.optimize import brentq
    

    Solving the IVP:

        def f(y,t):
            return 0.1920516*y*np.sin(y)
    
        t = np.linspace(0,12,100)
        y = odeint(f, 1, t).flatten()
        plt.plot(t,y);
    

    Finally, interpolating the solution and solving:

        t0 = brentq(lambda y: y_interp(y)-3, 0, 12)
        t0
    

    Resulting in:
    8.194529565299755

    mark
  • edited August 20

    @ChristianDonaldson , @andrew , @Austin , @StephenC

    I think you've all got some good work here but need a little help with formatting. Keep in mind that formatting is done with Markdown, as described in [this forum post]. In particular, a block of code should should be indented four spaces to appear as a code block in your post. For example, if we want to our post to look like:

    def f(x):
        a = 2
        return x+a
    f(3)
    
    # Output:
    # 5
    

    I'd type:

        def f(x):
            a = 2
            return x+a
        f(3)
    
        # Output:
        # 5
    

    Note that I've include the output as a Python comment.

    You folks are using backticks to indicate code but that is meant for inline code. For example, if I wanted to write: Use the def statement to define a function, I'd just wrap the "def" in backticks. Note that inline code doesn't preserve linebreaks, which i why you all are having problems with that.

  • @Dylan

    It looks like you're having a bit of difficulty with the mathematical typesetting. Please don't upload images as a replacement for the typesetting facilities.

  • edited August 20

    SAM: 19,1,13

    $r = .19113$
    $$y'=.19113ysin(y)$$

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    
    # A numerical rootfinder
    from scipy.optimize import brentq
    
    # Redefine f as a function of both y and t
    def f(y,t):
        return .19113*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,1000)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    odeint(f, 1, t)[:8]
    

    array([[1. ],
    [1.00322837],
    [1.00647386],
    [1.00973657],
    [1.01301659],
    [1.01631402],
    [1.01962897],
    [1.02296157]])

    odeint(f,1,t)[:8].flatten()
    

    array([1. , 1.00322837, 1.00647386, 1.00973657, 1.01301659,
    1.01631402, 1.01962897, 1.02296157])

    y_interp = interp1d(t,y)
    y_interp(20)
    

    array(3.1414658)

    y_interp(20)
    

    array(3.1414658)

    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    

    8.23309851814031 -This is our solution

    y_interp(t0)
    

    array(3.) -Checking

    plt.plot(t,y)
    plt.plot([t0,t0],[0,6],'m-')
    plt.plot([0,25],[3,3],'m-')
    plt.plot(t0,3,'co');
    

    mark
  • edited August 20

    Based on my forum username "Badler" my personal IVP is:
    $$y'=0.21412517ysin(y), y(0)=1$$
    Where $$ r=0.21412517 $$
    First I'm going to load the libraries:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    
    # A numerical rootfinder
    from scipy.optimize import brentq
    

    Now I'll set up the IVP:

    # Redefine f as a function of both y and t
    def f(y,t): 
        return 0.21412517*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,10,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    


    Now I will interpolate and solve:

    # y_interp(10)
    array(3.1168268)
    # brentq?
    t0 = brentq(lambda y: y_interp(y)-3, 0, 10)
    t0
    output=7.349488301836956
    
    mark
  • edited August 20

    My login name is Jules so my numbers are 10, 21, 12, 5, and 19. Thus, my $r=0.102112519.$ So, my ODE is $$ y'=0.102112519ysin(y).$$

    # A graphics library
    import matplotlib.pyplot as plt
    # Basic numerical tools
    import numpy as np
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    # A numerical rootfinder
    from scipy.optimize import brentq
    

    Now, let's set up and solve the IVP:

    Redefine f as a function of both y and t
    def f(y,t):
        return 0.102112519*y*np.sin(y)
    
    Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,100)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Finally, I'll interpolate the solution and solve it:

    y_interp = interp1d(t,y)
    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    
    #Output:
    15.411537406090577
    
    mark
  • maxmax
    edited August 20

    From "max" I obtained the IVP

    $$
    y' = 0.13124ysin(y), y(0)=1.
    $$

    First load the libraries:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d
    
    # A numerical rootfinder
    from scipy.optimize import brentq
    

    Then define the IVP and get a numeric solution:

    # Redefine f as a function of both y and t
    def f(y,t):
        return 0.13124*y*np.sin(y)
    
    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,20,1000)
    
    # Solve it!
    y = odeint(f, 1, t).flatten()
    
    # Plot it!!
    plt.plot(t,y);
    

    Now we interpolate the solution and then find $t$ in $y(t)=3$:

    y_interp = interp1d(t,y)
    t0 = brentq(lambda y: y_interp(y)-3, 0, 20)
    t0
    

    This code produces the output

    11.99016421201848

    So $y(t)=3$ at $t=11.99016421201848$

    This is further verified by running y_interp(t0), which outputs array(3.)

    mark
  • edited August 27

    My IVP is
    $$y^1= 0.191181ysin(y), y(0) = 1. $$

    First I need to load the libraries:

    import matplotlib.pyplot as plt
    import numpy as np
    from scipy.integrate import odeint
    from scipy.interpolate import interp1d
    from scipy.optimize import brentq
    

    Now I'm going to set up and solve the IVP:

    def f(y,t):
        return 0.191181*(y)*np.sin(y)
    t = np.linspace(0,20,100)
    y = odeint(f, 1, t).flatten()
    y
    plt.plot(t,y)
    

    My final output is found by the following code:

    y_interp = interp1d(t,y)
    y_interp
    t0 = brentq(lambda y: y_interp(y)-3,0,10)
    t0
    

    So our final output is 8.232968014338422

  • edited August 21

    My IVP is:

    $$y' = 0.1651818914y\sin y; y(0) = 1$$

    Load libraries:

    # A graphics library
    import matplotlib.pyplot as plt
    
    # Basic numerical tools
    import numpy as np
    
    # Scipy's numerical ODE solver
    from scipy.integrate import odeint
    
    # Tool to interpolate the solution
    from scipy.interpolate import interp1d
    
    # Tool to solve the interpolated solution for a specific value of t
    from scipy.optimize import brentq
    

    Set up and solve the IVP:

    # Redefine f as a function of both y and t
    def f(y,t):
        return 0.1651818914 * y * np.sin(y)
    t = np.linspace(0, 20, 100)
    y = odeint(f, 1, t).flatten()
    
    plt.plot(t,y)
    

    Interpolate the solution, solve, and visualize:

    y_interp = interp1d(t,y)
    t1 = brentq(lambda t: y_interp(t)-3,0,10)
    t1
    #Output: 9.527669205421008
    
    plt.plot(t,y)
    plt.plot([t1,t1],[0,5],'k--')
    plt.plot([0, 20], [3, 3], 'k--')
    plt.plot(t1, 3, 'ro')
    

    mark
  • Load Library

    # A graphics library
    import matplotlib.pyplot as plt

    # Basic numerical tools
    import numpy as np

    # Scipy's numerical ODE solver
    from scipy.integrate import odeint

    # Convert data generated by odeint to a computable function
    from scipy.interpolate import interp1d

    # A numerical rootfinder
    from scipy.optimize import brentq

    Set up and solve the equation

    # Redefine f as a function of both y and t
    def f(y,t):
    return 0.23977514152120ynp.sin(y)

    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,10,100)

    # Solve it!
    y = odeint(f, 1, t).flatten()

    # Plot it!!
    plt.plot(t,y);

    # Some t values that are appropriate to plot the solution
    t = np.linspace(0,10,100)

    # Solve it!
    y = odeint(f, 1, t).flatten()

    # Plot it!!
    plt.plot(t,y);

    Interpolate, solve and visualize

    y_interp = interp1d(t,y)
    y_interp
    t0 = brentq(lambda y: y_interp(y)-3, 0, 10)
    t0
    #output: 6.562862752442042
    plt.plot(t,y)
    plt.plot([t0,t0],[0,6],'k--')
    plt.plot([0,50],[3,3],'k--')
    plt.plot(t0,3,'ro');

  • @Wiggenout You should indent anything that you want to appear as a code block four spaces. Have a look at this post and the references therein.

Sign In or Register to comment.