Euler's method: General exploration
- Author:
- Juan Carlos Ponce Campuzano
- Topic:
- Differential Equation
1.3.1 to 1.3.2 Euler's Method for Solving ODEs Numerically, Euler's method using Matlab
Script
#-----------------------------------
# Define ODE y' = f
#-----------------------------------
f(t,y) = sin(1-t)*(y-1)*(y+1)
#-----------------------------------
# Define increment
#-----------------------------------
Δt = Slider(0.01, 1, 0.01, 1, 400, false, true, false, false)
SetValue(Δt, 0.05)
#-----------------------------------
# Euler's method
#-----------------------------------
n = 50
Execute( Join( {"T0 = 0"}, Zip( "T"+ (k+1) +" = T"+ k +" + Δt", k, 0..n ) ) )
Execute( Join( {"Y0 = 0.4"}, Zip( "Y"+ (k+1) +" = Y"+ k +"+ f(T"+ k +", Y"+ k +") * Δt", k, 0..n ) ) )
Execute( Join( {"P0 = (T0, Y0)"}, Zip( "P"+ (k+1) +" = (T"+ (k+1) +", Y"+ (k+1) +")", k, 0..n ) ) )
# 4th Degree Runge-Kutta solution with step = 0.001
sol = SolveODE(f, T0, Y0, 3, 0.001)
#-----------------------------------
# Settings
#-----------------------------------
Execute( Join( {"X0 = (T0, 0)"}, Zip( "X"+ (k+1) +" = (T"+ (k+1) +", 0)", k, 0..n ) ) )
# Texts
T = Zip("y_{"+k"}", k, 0..n)
U = Zip("t_{"+k"}", k, 0..n)
Execute(Zip("SetCaption(P"+ k +", Element(T, "+ (k+1) +"))", k, 0..n))
Execute(Zip("SetCaption(X"+ k +", Element(U, "+ (k+1) +"))", k, 0..n))
b = "Blue"
r = "Red"
Execute(Zip("SetColor(X"+ k +", b)", k, 0..n))
Execute(Zip("SetColor(P"+ k +", r)", k, 0..n))
Execute(Zip("SetPointSize(X"+ k +", 5)", k, 0..n))
Execute(Zip("SetPointSize(P"+ k +", 5)", k, 0..n))
Execute(Zip("SetPointStyle(X"+ k +", 10)", k, 0..n))
Execute(Zip("SetPointStyle(P"+ k +", 10)", k, 0..n))
Execute(Zip("S"+ k +" = Segment(P"+ k +", X"+ k +")", k, 0..n))
Execute(Zip("SetLineStyle(S"+ k +", 2)", k, 0..n))
Execute(Zip("SetLineThickness(S"+ k +", 2)", k, 0..n))
SL = Join({P0}, CellRange(P1, P50))
appx = Polyline(SL)
#-----------------------------------------------
Finally add in Update this script to Δt
#-----------------------------------------------
If( Δt <= 0.04, Execute(Zip("ShowLabel(X"+ k +", false)", k, 1..n)), Execute(Zip("ShowLabel(X"+ k +", true)", k, 1..n)))
If( Δt <= 0.04, Execute(Zip("ShowLabel(P"+ k +", false)", k, 1..n)), Execute(Zip("ShowLabel(P"+ k +", true)", k, 1..n)))