Solving differential equations numerically

Question

I tried solving a very simple equation f = t**2 numerically. I coded a for-loop, so as to use f for the first time step and then use the solution of every loop through as the inital function for the next loop.

I am not sure if my approach to solve it numerically is correct and for some reason my loop only works twice (one through the if- then the else-statement) and then just gives zeros.

Any help very much appreciatet. Thanks!!!

## IMPORT PACKAGES
import numpy as np
import math
import sympy as sym
import matplotlib.pyplot as plt

## Loop to solve numerically

for i in range(1,4,1):
    if i == 1:
        f_old = t**2
        print(f_old)
    else: 
        f_old = sym.diff(f_old, t).evalf(subs={t: i})
        f_new = f_old + dt * (-0.5 * f_old)
        f_old = f_new
        print(f_old)

Answer 1

Scipy.integrate package has a function called odeint that is used for solving differential equations

Here are some resources Link 1 Link 2

y = odeint(model, y0, t)

model: Function name that returns derivative values at requested y and t values as dydt = model(y,t)

y0: Initial conditions of the differential states

t: Time points at which the solution should be reported. Additional internal points are often calculated to maintain accuracy of the solution but are not reported.

Example that plots the results as well :

import numpy as np
from scipy.integrate import odeint
import matplotlib.pyplot as plt

# function that returns dy/dt
def model(y,t):
    k = 0.3
    dydt = -k * y
    return dydt

# initial condition
y0 = 5

# time points
t = np.linspace(0,20)

# solve ODE
y = odeint(model,y0,t)

# plot results
plt.plot(t,y)
plt.xlabel('time')
plt.ylabel('y(t)')
plt.show()