Python - Lambda Week: Solving Differential Equations with Runge Kutta 4th Order Method
Welcome to Python Week! This we we're going to cover calculus and the keyword lambda.
Note: All Python scripts presented this week were created using a TI-NSpire CX II CAS. As of June 2022, the lambda keyword is available on all calculators (in the United States) that have Python. If you are not sure, please check your calculator manual.
Solving Differential Equations
This following script solves the differential equation:
y' = dy/dx = f(x,y)
with initial condition y(x0) = y0
Repeat the steps for each step size h:
f1 = h * f(x0, y0)
f2 = h * f(x0 + h/2, y0 + f1/2)
f3 = h * f(x0 + h/2, y0 + f2/2)
f4 = h * f(x0 + h, y0 + f3)
x0 = x0 + h (update x)
y0 = y0 + (f1 + 2*f2 + 2*f3 + f4)/6 (update y)
The small h is, the more accurate the calculated coordinates are.
rk4lam.py: Runge Kutta 4th Order Method
All answers are stored in the nested list t.
from math import *
print("Runge Kutta 4th Order")
print("Math Module imported")
f=eval("lambda x,y:"+input("dy/dx = "))
# must call for float numbers one at a time
x0=eval(input("x0 = "))
y0=eval(input("y0 = "))
h=eval(input("h = "))
# ask for an integer
n=int(input("number of steps: "))
# set up table
t=[[x0,y0]]
# main loop
for i in range(n):
f1=h*f(x0,y0)
f2=h*f(x0+h/2,y0+f1/2)
f3=h*f(x0+h/2,y0+f2/2)
f4=h*f(x0+h,y0+f3)
x0=x0+h
y0=y0+(f1+2*f2+2*f3+f4)/6
print([x0,y0])
t.append([x0,y0])
print("Done. Recall t for table.")
Examples
Results are rounded to five digits.
Example 1:
dy/dx = 2*x*y + x, y(0) = 0, h = 0.1, 5 steps
(Real solution: y = 1/2 * (e^(x^2) - 1))
Results (which matches the exact results):
x = 0.1, y ≈ 0.00503
x = 0.2, y ≈ 0.02041
x = 0.3, y ≈ 0.04709
x = 0.4, y ≈ 0.08676
x = 0.5, y ≈ 0.14201
Example 2:
dy/dx = ln x + y, y(10) = 1
(Real Solution: y = [∫(ln t * e^(-t) dt, t = 10 to x) + e^(-10)] * e^x
Exact Results:
x = 11, y ≈ 6.74551
x = 12, y ≈ 22.51732
x = 13, y ≈ 65.53659
x = 14, y ≈ 182.60824
x = 15, y ≈ 500.96552
Runge Kutta with h = 1, 5 steps:
x = 11, y ≈ 6.71066
x = 12, y ≈ 22.33376
x = 13, y ≈ 64.78988
x = 14, y ≈ 179.90761
x = 15, y ≈ 491.80768
Runge Kutta with h = 0.1, 50 steps:
x = 11, y ≈ 6.74551 (recall t[10])
x = 12, y ≈ 22.51728 (t[20])
x = 13, y ≈ 65.53643 (t[30])
x = 14, y ≈ 182.60766 (t[40])
x = 15, y ≈ 500.96358 (t[50])
This ends Python week for now, I hope you find this week helpful and resourceful.
Until next time,
Eddie
All original content copyright, © 2011-2022. Edward Shore. Unauthorized use and/or unauthorized distribution for commercial purposes without express and written permission from the author is strictly prohibited. This blog entry may be distributed for noncommercial purposes, provided that full credit is given to the author.
