Added demonstration codes for IVP methods

drreynolds · Apr 5, 2021 · 9a8639b · 9a8639b
1 parent a5e981c
commit 9a8639b
Show file tree

Hide file tree

Showing 13 changed files with 1,053 additions and 0 deletions.
diff --git a/InitialValue/AB2_step.m b/InitialValue/AB2_step.m
@@ -0,0 +1,27 @@
+function [unew,fcur] = AB2_step(f, tcur, ucur, h, fold)
+% usage: [unew,fcur] = AB2_step(f, tcur, ucur, h, fold)
+%
+% Adams-Bashforth method of order 2 for one step of the ODE
+% problem,
+%    u' = f(t,u), t in tspan,
+%    u(t0) = u0.
+%
+% Inputs:  f = function for ODE right-hand side, f(t,u)
+%          tcur = current time
+%          ucur = current solution
+%          h = time step size
+%          fold = RHS evaluated at previous time step
+% Outputs: unew = updated solution
+%          fcur = RHS evaluated at current time step
+%
+% Daniel R. Reynolds
+% SMU Mathematics
+% Math 4315
+
+% get ODE RHS at this time step
+fcur = f(tcur, ucur);
+
+% update solution in time
+unew = ucur + h/2*(3*fcur - fold);
+
+% end of function
diff --git a/InitialValue/Heun_step.m b/InitialValue/Heun_step.m
@@ -0,0 +1,32 @@
+function unew = Heun_step(f, told, uold, h)
+% usage: unew = Heun_step(f, told, uold, h)
+%
+% Heun's Runge-Kutta of order 2 solver for one step of the
+% ODE problem,
+%    u' = f(t,u), t in tspan,
+%    u(t0) = u0.
+%
+% Inputs:  f = function for ODE right-hand side, f(t,u)
+%          told = current time
+%          uold = current solution
+%          h = time step size
+% Outputs: unew = updated solution
+%
+% Daniel R. Reynolds
+% SMU Mathematics
+% Math 4315
+
+% call f to get ODE RHS at this time step
+f1 = f(told, uold);
+
+% set z1 and z2
+z1 = uold;
+z2 = uold + h*f1;
+
+% get f(t+h,z2)
+f2 = f(told+h, z2);
+
+% update solution in time
+unew = uold + h/2*(f1+f2);
+
+% end of function
diff --git a/InitialValue/RK2_step.m b/InitialValue/RK2_step.m
@@ -0,0 +1,32 @@
+function unew = RK2_step(f, told, uold, h)
+% usage: unew = RK2_step(f, told, uold, h)
+%
+% Runge-Kutta of order 2 solver for one step of the ODE
+% problem,
+%    u' = f(t,u), t in tspan,
+%    u(t0) = u0.
+%
+% Inputs:  f = function for ODE right-hand side, f(t,u)
+%          told = current time
+%          uold = current solution
+%          h = time step size
+% Outputs: unew = updated solution
+%
+% Daniel R. Reynolds
+% SMU Mathematics
+% Math 4315
+
+% get ODE RHS at this time step
+f1 = f(told, uold);
+
+% set z1 and z2
+z1 = uold;
+z2 = uold + h/2*f1;
+
+% get f(t+h/2,z2)
+f2 = f(told+h/2, z2);
+
+% update solution in time
+unew = uold + h*f2;
+
+% end of function
diff --git a/InitialValue/bwd_Euler_step.m b/InitialValue/bwd_Euler_step.m
@@ -0,0 +1,33 @@
+function unew = bwd_Euler_step(f, f_u, told, uold, h)
+% Usage: unew = bwd_Euler_step(f, f_u, told, uold, h)
+%
+% Backward Euler solver for one step of the ODE problem,
+%    u' = f(t,u), t in tspan,
+%    u(t0) = u0.
+%
+% Inputs:  f = function for ODE right-hand side, f(t,u)
+%          f_u = function for ODE Jacobian, f_u(t,u)
+%          told = current time
+%          uold = current solution
+%          h = time step size
+% Outputs: unew = updated solution
+%
+% Daniel R. Reynolds
+% SMU Mathematics
+% Math 4315
+
+% set nonlinear solver parameters
+maxit = 20;
+rtol = 1e-9;
+atol = 1e-12;
+output = false;
+
+% create implicit residual and Jacobian functions
+I = eye(length(uold));
+F = @(unew) unew - uold - h*f(told+h, unew);
+A = @(unew) I - h*f_u(told+h, unew);
+
+% perform implicit solve
+[unew, its] = newton(F, A, uold, maxit, rtol, atol, output);
+
+% end of function
diff --git a/InitialValue/driver1.m b/InitialValue/driver1.m
@@ -0,0 +1,131 @@
+% Driver to test various ODE solvers for the test problem
+%    u' = (u+t^2-2)/(t+1),   t in [0,5]
+%    u(0) = 2.
+%
+% Daniel R. Reynolds
+% SMU Mathematics
+% MATH 4315
+clear
+
+% ODE RHS and derivatives, initial condition, etc.
+f   = @(t,u,p) (u+t.^2-2)./(t+1);
+f_u = @(t,u,p) 1/(t+1);
+utrue = @(t) t.^2 + 2*t + 2 - 2*(t+1).*log(t+1);
+u0  = 2;
+t0  = 0;
+Tf  = 5;
+hvals = 0.5.^(2:8);
+
+% initialize error storage arrays
+err_fE     = zeros(size(hvals));
+err_RK2    = zeros(size(hvals));
+err_Heun   = zeros(size(hvals));
+err_AB2    = zeros(size(hvals));
+err_bE     = zeros(size(hvals));
+err_trap   = zeros(size(hvals));
+
+% iterate over our h values
+for j=1:length(hvals)
+
+   % create time output array
+   h = hvals(j);
+   N = floor((Tf-t0)/h);
+   tvals = linspace(t0,Tf,N+1);
+
+   % initialize solutions
+   u_fE     = u0;
+   u_RK2    = u0;
+   u_Heun   = u0;
+   u_AB2    = u0;
+   u_bE     = u0;
+   u_trap   = u0;
+
+   % try out our methods
+   for i=1:N
+
+      t = tvals(i);
+      tnew = t + h;
+
+      u_fE = fwd_Euler_step(f, t, u_fE, h);
+      err_fE(j) = max([err_fE(j), abs(u_fE-utrue(tnew))]);
+
+      u_RK2 = RK2_step(f, t, u_RK2, h);
+      err_RK2(j) = max([err_RK2(j), abs(u_RK2-utrue(tnew))]);
+
+      u_Heun = Heun_step(f, t, u_Heun, h);
+      err_Heun(j) = max([err_Heun(j), abs(u_Heun-utrue(tnew))]);
+
+      % for first AB2 step, store fold and use RK2 for step; subsequently use AB2
+      if (i == 1)
+         fold_AB2 = f(t, u_AB2);
+         u_AB2 = RK2_step(f, t, u_AB2, h);
+      else
+         [u_AB2,fold_AB2] = AB2_step(f, t, u_AB2, h, fold_AB2);
+      end
+      err_AB2(j) = max([err_AB2(j), abs(u_AB2-utrue(tnew))]);
+
+      u_bE = bwd_Euler_step(f, f_u, t, u_bE, h);
+      err_bE(j) = max([err_bE(j), abs(u_bE-utrue(tnew))]);
+
+      u_trap = trap_step(f, f_u, t, u_trap, h);
+      err_trap(j) = max([err_trap(j), abs(u_trap-utrue(tnew))]);
+
+   end
+end
+
+% output convergence results
+disp('Results for Fwd Euler:')
+err = err_fE;
+fprintf('   h = %10g,  err = %.2e\n',hvals(1),err(1))
+for i=2:length(hvals)
+   fprintf('   h = %10g,  err = %.2e,  rate = %g\n',...
+           hvals(i),err(i),log(err(i)/err(i-1))/log(hvals(i)/hvals(i-1)))
+end
+
+disp('Results for RK2:')
+err = err_RK2;
+fprintf('   h = %10g,  err = %.2e\n',hvals(1),err(1))
+for i=2:length(hvals)
+   fprintf('   h = %10g,  err = %.2e,  rate = %g\n',...
+           hvals(i),err(i),log(err(i)/err(i-1))/log(hvals(i)/hvals(i-1)))
+end
+
+disp('Results for Heun:')
+err = err_Heun;
+fprintf('   h = %10g,  err = %.2e\n',hvals(1),err(1))
+for i=2:length(hvals)
+   fprintf('   h = %10g,  err = %.2e,  rate = %g\n',...
+           hvals(i),err(i),log(err(i)/err(i-1))/log(hvals(i)/hvals(i-1)))
+end
+
+disp('Results for AB2:')
+err = err_AB2;
+fprintf('   h = %10g,  err = %.2e\n',hvals(1),err(1))
+for i=2:length(hvals)
+   fprintf('   h = %10g,  err = %.2e,  rate = %g\n',...
+           hvals(i),err(i),log(err(i)/err(i-1))/log(hvals(i)/hvals(i-1)))
+end
+
+disp('Results for Bwd Euler:')
+err = err_bE;
+fprintf('   h = %10g,  err = %.2e\n',hvals(1),err(1))
+for i=2:length(hvals)
+   fprintf('   h = %10g,  err = %.2e,  rate = %g\n',...
+           hvals(i),err(i),log(err(i)/err(i-1))/log(hvals(i)/hvals(i-1)))
+end
+
+disp('Results for Trapezoidal:')
+err = err_trap;
+fprintf('   h = %10g,  err = %.2e\n',hvals(1),err(1))
+for i=2:length(hvals)
+   fprintf('   h = %10g,  err = %.2e,  rate = %g\n',...
+           hvals(i),err(i),log(err(i)/err(i-1))/log(hvals(i)/hvals(i-1)))
+end
+
+% display true solution
+figure()
+plot(tvals, utrue(tvals))
+xlabel('t')
+ylabel('u(t)')
+
+% end of script
diff --git a/InitialValue/driver1.py b/InitialValue/driver1.py
@@ -0,0 +1,131 @@
+#!/usr/bin/env python3
+#
+# Script to test various ODE solvers for the test problem
+#    u' = (u+t^2-2)/(t+1),   t in [0,5]
+#    u(0) = 2.
+#
+# Daniel R. Reynolds
+# SMU Mathematics
+# Math 4315
+
+# imports
+import numpy
+import matplotlib.pyplot as plt
+from ivp_methods import *
+
+# ODE RHS and Jacobian, initial condition, etc.
+def f(t,u):
+    return (u+t*t-2)/(t+1)
+def f_u(t,u):
+    return (1.0/(t+1))
+def utrue(t):
+    return (t**2 + 2*t + 2 - 2*(t+1)*numpy.log(t+1))
+u0  = 2.0
+t0  = 0.0
+Tf  = 5.0
+hvals = [0.5**2, 0.5**3, 0.5**4, 0.5**5, 0.5**6, 0.5**7, 0.5**8]
+
+# initialize error storage arrays
+err_fE     = numpy.zeros((len(hvals),1), dtype=float)
+err_RK2    = numpy.zeros((len(hvals),1), dtype=float)
+err_Heun   = numpy.zeros((len(hvals),1), dtype=float)
+err_AB2    = numpy.zeros((len(hvals),1), dtype=float)
+err_bE     = numpy.zeros((len(hvals),1), dtype=float)
+err_trap   = numpy.zeros((len(hvals),1), dtype=float)
+
+# iterate over our h values
+for j, h in enumerate(hvals):
+
+    # create time output array
+    N = int((Tf-t0)/h)
+    tvals = numpy.linspace(t0, Tf, N+1)
+
+    # initialize solutions
+    u_fE     = u0
+    u_RK2    = u0
+    u_Heun   = u0
+    u_AB2    = u0
+    u_bE     = u0
+    u_trap   = u0
+
+    # try out our methods
+    for i in range(N):
+
+        t = tvals[i]
+        tnew = t + h
+
+        u_fE = fwd_Euler_step(f, t, u_fE,  h)
+        err_fE[j] = max(err_fE[j], abs(u_fE-utrue(tnew)))
+
+        u_RK2 = RK2_step(f, t, u_RK2, h)
+        err_RK2[j] = max(err_RK2[j], abs(u_RK2-utrue(tnew)))
+
+        u_Heun = Heun_step(f, t, u_Heun, h)
+        err_Heun[j] = max(err_Heun[j], abs(u_Heun-utrue(tnew)))
+
+        # for first AB2 step, store fold and use RK2 for step; subsequently use AB2
+        if (i == 0):
+            fold_AB2 = f(t, u_AB2)
+            u_AB2 = RK2_step(f, t, u_AB2, h)
+        else:
+            u_AB2, fold_AB2 = AB2_step(f, t, u_AB2, h, fold_AB2)
+        err_AB2[j] = max(err_AB2[j], abs(u_AB2-utrue(tnew)))
+
+        u_bE = bwd_Euler_step(f, f_u, t, u_bE,  h)
+        err_bE[j] = max(err_bE[j], abs(u_bE-utrue(tnew)))
+
+        u_trap = trap_step(f, f_u, t, u_trap,  h)
+        err_trap[j] = max(err_trap[j], abs(u_trap-utrue(tnew)))
+
+
+# output convergence results
+print('Results for Fwd Euler:')
+err = err_fE
+print("   h = %10g,  err = %.2e" % (hvals[0], err[0]))
+for i in range(1,len(hvals)):
+    print("   h = %10g,  err = %.2e,  rate = %g" %
+          (hvals[i], err[i], numpy.log(err[i]/err[i-1])/numpy.log(hvals[i]/hvals[i-1])) )
+
+print('Results for RK2:')
+err = err_RK2
+print("   h = %10g,  err = %.2e" % (hvals[0], err[0]))
+for i in range(1,len(hvals)):
+    print("   h = %10g,  err = %.2e,  rate = %g" %
+          (hvals[i], err[i], numpy.log(err[i]/err[i-1])/numpy.log(hvals[i]/hvals[i-1])) )
+
+print('Results for Heun:')
+err = err_Heun
+print("   h = %10g,  err = %.2e" % (hvals[0], err[0]))
+for i in range(1,len(hvals)):
+    print("   h = %10g,  err = %.2e,  rate = %g" %
+          (hvals[i], err[i], numpy.log(err[i]/err[i-1])/numpy.log(hvals[i]/hvals[i-1])) )
+
+print('Results for AB2:')
+err = err_AB2
+print("   h = %10g,  err = %.2e" % (hvals[0], err[0]))
+for i in range(1,len(hvals)):
+    print("   h = %10g,  err = %.2e,  rate = %g" %
+          (hvals[i], err[i], numpy.log(err[i]/err[i-1])/numpy.log(hvals[i]/hvals[i-1])) )
+
+print('Results for Bwd Euler:')
+err = err_bE
+print("   h = %10g,  err = %.2e" % (hvals[0], err[0]))
+for i in range(1,len(hvals)):
+    print("   h = %10g,  err = %.2e,  rate = %g" %
+          (hvals[i], err[i], numpy.log(err[i]/err[i-1])/numpy.log(hvals[i]/hvals[i-1])) )
+
+print('Results for Trapezoidal:')
+err = err_trap
+print("   h = %10g,  err = %.2e" % (hvals[0], err[0]))
+for i in range(1,len(hvals)):
+    print("   h = %10g,  err = %.2e,  rate = %g" %
+          (hvals[i], err[i], numpy.log(err[i]/err[i-1])/numpy.log(hvals[i]/hvals[i-1])) )
+
+# display true solution
+plt.figure()
+plt.plot(tvals, utrue(tvals))
+plt.xlabel('$t$')
+plt.ylabel('$u(t)$')
+plt.show()
+
+# end of script