Numerical and Graphical Solutions of Volterra Integral Equations of the Second KindLeigh C. Becker and Micah WheelerDepartment of Mathematics Christian Brothers University 650 E. Parkway South Memphis, TN 38104 lbecker@cbu.edu

The principal part of this worksheet is the Maple procedure in Section 4 called Volterra. It can be used to numerically solve a Volterra integral equation of the second kind, namely, an equation of the form (1.1) NiMvLUkieEc2IjYjSSJ0R0YmLCYtSSJmR0YmRiciIiItSSRJbnRHNiRJKnByb3RlY3RlZEdGMEkoX3N5c2xpYkdGJjYkLUkia0dGJjYlRihJInNHRiYtRiU2I0Y2L0Y2OyIiIUYoRiw=, where x, f, and k are vector-valued functions with m components. If f and k are continuous and k satisfies a Lipschitz condition with respect to x, then a unique solution x(t) of (1.1) exists. For more information, see Burton [2, secs. 2.1 and 3.2] or Linz [4, sec. 4.1]. Since (1.1) is a vector equation, it is equivalent to a system consisting of the m scalar equations (1.2) NiMvLSZJInhHNiI2I0kiaUdGJzYjSSJ0R0YnLCYtJkkiZkdGJ0YoRioiIiItSSRJbnRHRic2JC0mSSJrR0YnRig2JUYrSSJzR0YnLUYmNiNGOC9GODsiIiFGK0Yw (i = 1, ... , m). Volterra approximates the solution x(t) at equally spaced points tn using the implicit trapezoidal rule in conjunction with Newton's method for nonlinear systems. Computer programs in other programming languages that will numerically solve (1.1) using this approach already exist in the literature. Among them is a Pascal program in [4, ch. 12]. There is a Fortran program in [5, ch. 18.2] that will numerically solve the linear version of (1.1), i.e., when k(t, s, x(s)) = k(t, s) x(s), with the explicit trapezoidal rule. Despite the availability of these and other such programs, there are advantages in using the Volterra procedure; notably among them are: The tedium of computing all of the first-order partial derivatives of the components of the kernel k by hand (cf. equation (2.9) in Section 2) is eliminated with the Jacobian and unapply commands. Newton's method for nonlinear systems entails solving linear algebraic systems (cf. (2.11) in Section 2). This is achieved with the command LinearSolve.With the graphics facilities of Maple, the numerical solutions obtained with Volterra are easily graphed with the pointplot or pointplot3d commands. Four examples using Volterra to numerically solve equations of the form (1.1)--or those that can be written in that form--and graph the results are included in this worksheet. The first example is a system of two nonlinear equations taken from Linz [4, p. 201]; it is considered in Section 5. In Section 6, we use Volterra to numerically solve the scalar linear Volterra equation considered by Jerri [3, pp. 90-93]. In Section 7, we choose a scalar Volterra integro-differential equation of the form (1.3) NiMvLUklZGlmZkdJKnByb3RlY3RlZEdGJjYkLUkieEc2IjYjSSJ0R0YqRiwsJiomLUkiYUdGKkYrIiIiRihGMUYxLUkkSW50R0YqNiQqJi1JImJHRio2JEYsSSJzR0YqRjEtRik2I0Y5RjEvRjk7IiIhRixGMQ== from Example 1.2.1 in Burton [2, p. 8] to show how to change such an equation to an equivalent system of Volterra integral equations of the second kind so that we can numerically solve it with Volterra. The last example is found in Section 8, where we concocted a system of three nonlinear Volterra integral equations.

The Volterra procedure approximates the solution x(t) of (1.1) at the equally spaced points tn = t0 + nh for n = 1, ... , N, where t0 = 0 and N is the total number of steps of size h. Xn denotes the approximation of x(t) at t = tn. In the Volterra procedure itself, tn and Xn are denoted by t[n] and X[n], respectively. Let us briefly describe the mathematics underlying the procedure. Setting t = tn in (1.1), we have(2.1) NiMvLUkieEc2IjYjJkkidEdGJjYjSSJuR0YmLCYtSSJmR0YmRiciIiItSSRpbnRHRiY2JC1JImtHRiY2JUYoRiktRiU2I0YpL0YpOyIiIUYoRi8=.By the composite trapezoidal rule (cf. [1, p. 200]) , an approximation of the integral in (2.1) is (2.2) NiMqJkkiaEc2IiIiIiIiIyEiIg== [NiMtSSJrRzYiNiUmSSJ0R0YlNiNJIm5HRiUmRig2IyIiIS1JInhHRiU2I0Yr + 2 NiMtSSRTdW1HNiI2JC1JImtHRiU2JSZJInRHRiU2I0kibkdGJSZGKzYjSSJqR0YlLUkieEdGJTYjRi4vRjA7IiIiLCZGLUY2RjYhIiI= + NiMtSSJrRzYiNiUmSSJ0R0YlNiNJIm5HRiVGJy1JInhHRiU2I0Yn].Replacing x(tn) in (2.1) and (2.2) by Xn, we obtain the implicit trapezoidal rule (2.3) NiMvJkkiWEc2IjYjSSJuR0YmLUkiZkdGJjYjJkkidEdGJkYn + NiNJImhHNiI=[ NiMqJiIiIkYkIiIjISIi NiMtSSJrRzYiNiUmSSJ0R0YlNiNJIm5HRiUmRig2IyIiISZJIlhHRiVGLA== + NiMtSSRTdW1HNiI2JC1JImtHRiU2JSZJInRHRiU2I0kibkdGJSZGKzYjSSJqR0YlJkkiWEdGJUYvL0YwOyIiIiwmRi1GNUY1ISIi + NiMqJiIiIkYkIiIjISIi NiMtSSJrRzYiNiUmSSJ0R0YlNiNJIm5HRiVGJyZJIlhHRiVGKQ==],where X0 = f(0) since x(0) = f(0). Defining NiMmSSZzaWdtYUc2IjYjSSJuR0Yl by(2.4) NiMmSSZzaWdtYUc2IjYjSSJuR0Yl := NiMtSSJmRzYiNiMmSSJ0R0YlNiNJIm5HRiU= + NiNJImhHNiI=[NiMqJiIiIkYkIiIjISIi NiMtSSJrRzYiNiUmSSJ0R0YlNiNJIm5HRiUmRig2IyIiISZJIlhHRiVGLA== + NiMtSSRTdW1HNiI2JC1JImtHRiU2JSZJInRHRiU2I0kibkdGJSZGKzYjSSJqR0YlJkkiWEdGJUYvL0YwOyIiIiwmRi1GNUY1ISIi ], we can rewrite (2.3) as(2.5) NiMmSSJYRzYiNiNJIm5HRiU= - NiMqJiIiIkYkIiIjISIi NiNJImhHNiI=NiMvLCYtSSJrRzYiNiUmSSJ0R0YnNiNJIm5HRidGKSZJIlhHRidGKyIiIiZJJnNpZ21hR0YnRishIiIiIiE=,where 0 denotes the zero vector. From (2.5), we see that Xn is the solution of the vector equation(2.6) NiMvLUkkcGhpRzYiNiNJInVHRiYiIiE=,where NiNJJHBoaUc2Ig== is the vector-valued function(2.7) NiMtSSRwaGlHNiI2I0kidUdGJQ== := NiNJInVHNiI= - NiMqJiIiIkYkIiIjISIi NiMsJiomSSJoRzYiIiIiLUkia0dGJjYlJkkidEdGJjYjSSJuR0YmRitJInVHRiZGJ0YnJkkmc2lnbWFHRiZGLSEiIg==. We will obtain an approximation to the solution Xn of (2.6) by way of the matrix-valued function G defined in (2.8). If A(u) is an m by m matrix-valued function that is invertible in a neighborhood of Xn, then Xn is a fixed point of (2.8) NiM+LUkiR0c2IjYjSSJ1R0YmLCZGKCIiIiomLUkiQUdGJkYnRiotSSRwaGlHRiZGJ0YqISIi.Assuming the components of G(u) have continuous first and second order partial derivatives and that the first order partial derivatives at Xn are equal to zero, it can be shown that if A(u) is set equal to the Jacobian matrix of the function NiNJJHBoaUc2Ig==, the iterates Xn(p) defined by (2.10) below will usually converge quadratically to Xn, provided the starting value is sufficiently close to Xn. The Jacobian matrix of NiNJJHBoaUc2Ig== is the m by m matrix J(u) with the element (2.9) J(u)ij = NiMtSSVEaWZmRzYiNiQtJkkkcGhpR0YlNiNJImlHRiU2I0kidUdGJSZGLTYjSSJqR0Yl = NiNJJmRlbHRhRzYiij - NiMqJiIiIkYkIiIjISIi NiNJImhHNiI= NiMtSSVEaWZmRzYiNiQtJkkia0dGJTYjSSJpR0YlNiUmSSJ0R0YlNiNJIm5HRiVGLUkidUdGJSZGMTYjSSJqR0Yl in row i and column j, where NiNJJmRlbHRhRzYiij is the Kronecker delta. Details of the statements made here follow from the discussion of Newton's method for nonlinear systems in [1, sec. 10.2]. Linz gives a brief outline of the trapezoidal rule and Newton's method for Volterra integral systems of the second kind in Section 12.1 of [4]. We obtain Xn from Xn-1 by setting Xn(0) := Xn-1 and then generating the iterates Xn(p) from (2.10) Xn(p) = G(Xn(p-1)) = Xn(p-1) - J -1(Xn(p-1))NiNJJHBoaUc2Ig==(Xn(p-1))for p = 1, 2, 3, ... . (This is Newton's method for nonlinear systems.) Let y denote the solution of the matrix equation(2.11) J(Xn(p-1))y = NiNJJHBoaUc2Ig==(Xn(p-1)). Then the iteration formula (2.10) becomes(2.12) Xn(p) = Xn(p-1) - y.We compute the solution y = J -1(Xn(p-1))NiNJJHBoaUc2Ig==(Xn(p-1)) using the command LinearSolve. The iterates Xn(p) are computed until the infinity norm of the vector y is less than a prescribed tolerance Tol. Then Xn is assigned the value of the last iterate. However, should the number of iterations at any stage exceed some assigned maximum value that we call Iter_Max, the procedure will be terminated.

restart: with(LinearAlgebra, LinearSolve, VectorNorm); with (plots, pointplot, display, pointplot3d, display3d); with(VectorCalculus, Jacobian); Digits := 10;

To implement Volterra, assign values to m, h, N, Tol, Iter_Max, respectively, wherem := number of equations; also, the number of unknown scalar functions xi(t);h := step size;N := total number of steps of size h;Tol := tolerance;Iter_Max := maximum number of Newton iterations. Volterra := proc(m::posint, h::positive, N::posint, Tol::positive, Iter_Max::posint) local i, j, p, Matrix_Partials; global n, t, k, kpartial, X, sigma, J, phi, y, epsilon; Let y, NiNJJnNpZ21hRzYi, NiNJJHBoaUc2Ig==, and Xn be vectors with m components: y := Vector(m); sigma := Vector(m); phi := Vector(m); for n from 0 to N do X[n] := Vector(m) NiI= end do; J is the m by m Jacobian matrix defined by equation (2.9): J := Matrix(m, m); The matrix kpartial consists of the first-order partial derivatives of k(t, s, x): Matrix_Partials := Jacobian(k(t, s, v), vlist); kpartial := unapply(Matrix_Partials, t, s, x); t[0] := 0; (t0 = 0 and X0 = f(0)) X[0] := f(t[0]); for n from 1 to N do t[n] := t[n-1] + h; sigma := f(t[n]) + 0.5*h*k(t[n], t[0], X[0]); (cf. (2.4)) for j from 1 to n-1 do sigma := sigma + h*k(t[n], t[j], X[j]) end do; X[n] := X[n-1]; p := 0; epsilon := 1; (To begin, choose a value for epsilon greater than Tol.) while (epsilon >= Tol) and (p <= Iter_Max) do p := p + 1; for i from 1 to m do (cf. (2.9)) for j from 1 to m do if i = j then J[i, j] := 1 - 0.5*h*kpartial(t[n], t[n], X[n])[i, j] else J[i, j] := - 0.5*h*kpartial(t[n], t[n], X[n])[i, j] end if; end do; end do; phi := X[n] - 0.5*h*k(t[n], t[n], X[n]) - sigma; (cf. (2.7)) LinearSolve(J, phi) computes the solution y of the linear system Jy = NiNJJHBoaUc2Ig==: y := LinearSolve(J, phi); VectorNorm(y) computes the infinity norm of y (cf. Sec. 2 and [1, p.419]). epsilon := VectorNorm(y); X[n] := X[n] - y; (cf. (2.12)) end do; (end of the "while" loop) Test for convergence: if epsilon >= Tol then "The maximum number of Newton iterations has been exceeded."; break end if; end do; (end of the "for n from 1 to N do" loop) end proc:

Linz [4, ch. 13] computes the numerical solution of the Volterra system NiMtJkkieEc2IjYjIiIiNiNJInRHRiY= = NiNJInRHNiI= +NiMtSSRpbnRHNiI2JDwjLCYqJiIiIkYqLCZGKkYqKiQtJkkieEdGJTYjRio2I0kic0dGJSIiI0YqISIiRioqJi1JJGNvc0dGJTYjLCZJInRHRiVGKkYyRjRGKi0mRi82I0YzRjFGKkYqL0YyOyIiIUY6 (5.1) NiMvLSZJInhHNiI2IyIiIzYjSSJ0R0YnLCYiIiJGLS1JJEludEdGJzYkPCMsJiomLUkkc2luR0YnNiMsJkYrRi1JInNHRichIiJGLS0mRiY2I0YtNiNGOEYtRi0qJi1JJGV4cEdGJzYjLCZGOEYtRitGOUYtLUYlRj1GLUY5L0Y4OyIiIUYrRi0= using the step size h = 0.2. See equations (13.1) and (13.2) on p. 201 in [4]. We use Volterra to numerically solve and graph (5.1), after which we compare the results with the values listed in Table 13.1 on p. 202 in Linz.Since system (5.1) consists of two equations, m = 2:m:=2; Enter the function f(t):f := proc(t) Vector([t, 1]) end proc: Check that f has been entered correctly:f(t); Enter the kernel k(t, s, x): k := proc(t, s, x::Vector) Vector([1/(1 + x[1]^2) + cos(t-s)*x[2], sin(t-s)*x[1] - exp(s-t)*x[2]]) end proc: Check that k(t, s, x) has been entered correctly:v := Vector(m, symbol = x): vlist := convert(v, 'list'): k(t, s, v);Check that (5.1) has been entered correctly:for i from 1 to m do x[i](t) = f(t)[i] + Int(k(t, s, v)[i], s = 0..t) end do;Applying Volterra to (5.1) with step size h = 0.2:Volterra(2, 0.2, 50, 10^(-8),5);N:=n-1; (total number of steps of size h)d:=1; (In the following table, every value is displayed if d = 1; every other value if d = 2, etc.) M:=5; (M + 1 is the number of values that will be displayed in each column.) Numerical solution: Remark. In the following table, x1(t) and x2(t) denote the respective approximations of x1(t) and x2(t) for h = 0.2. (The approximations agree with the values Linz [4] gives on p. 202.)printf(" t x1(t) x2(t) \n"); for j from 0 to M do printf("%10.1f %10.5f %10.5f\n", t[d*j], X[d*j][1], X[d*j][2]) end do; Graphs:x1(t) vs. t: P[1] := pointplot([seq([t[n], X[n][1]], n = 0..N)], color=black, labels=[t, "x1"], labelfont=[times, bold, 12]): P[2] := pointplot([seq([t[n], X[n][1]], n = 0..N)], color=magenta, style=line, thickness=1, labels=[t, "x1"], labelfont=[times, bold, 12]): display(P[1], P[2]); x2(t) vs. t:P[3] := pointplot([seq([t[n], X[n][2]], n = 0..N)], color=black, labels=[t, "x2"], labelfont=[times, bold, 12]): P[4] := pointplot([seq([t[n], X[n][2]], n = 0..N)], color=green, style=line, thickness=1, labels=[t, "x2"], labelfont=[times, bold, 12]): display(P[3], P[4]); x1(t) and x2(t) vs. t:P[5] := pointplot3d([seq([t[n], X[n][1], X[n][2]], n = 0..N)], view = [0..10, 0..10, 0...10], color=black, labels=[t, "x1", "x2"], labelfont=[times, bold, 12], axes = normal, tickmarks = [3,3,3], orientation = [45, 100]): P[6] := pointplot3d([seq([t[n], X[n][1], X[n][2]], n = 0..N)], view = [0..10, 0..10, 0...10], color=blue, labels=[t, "x1", "x2"], labelfont=[times, bold, 12], axes = normal, style=line, tickmarks = [3, 3, 3], thickness=1): display(P[5], P[6]);

We numerically solve the scalar Volterra equation(6.1) NiMvLUkieEc2IjYjSSJ0R0YmLCZGKCIiIi1JJEludEc2JEkqcHJvdGVjdGVkR0YuSShfc3lzbGliR0YmNiQqJiwmRihGKkkic0dGJiEiIkYqLUYlNiNGM0YqL0YzOyIiIUYoRjQ= and compare the results to the analytical solution x(t) = sin t. See Jerri [3, p. 78; pp. 90-93]. For the scalar equation (6.1), m = 1:m := 1;Enter the function f(t):f := proc(t) Vector([t]) end proc:Check that f has been entered correctly :f(t);Enter the kernel k(t, s, x):k := proc(t, s, x::Vector) Vector([(s - t)*x[1]]) end proc:Check that k(t, s, x) has been entered correctly:v := Vector(m, symbol = x): vlist := convert(v, 'list'): k(t, s, v);Check that (6.1) has been correctly entered:for i from 1 to m do x[i](t) = f(t)[i] + Int(k(t, s, v)[i], s = 0..t) end do; Applying Volterra to (6.1) with step size h = 1:Volterra(1, 1, 4, 10^(-5), 5);N:=n-1; (total number of steps of size h)d:=1; (In the following table, every value is displayed if d = 1; every other value if d = 2, etc.) M:=4; (M + 1 is the number of values that will be displayed in each column.) Numerical solution: Remark. The following approximations x1(t) of x1(t) = x(t) = sin t agree with the values Jerri [3] gives in Table 3.1 on p. 91.printf(" t x1(t) sin(t) \n"); for j from 0 to M do printf("%10.1f %10.4f %10.4f\n", t[d*j], X[d*j][1], sin(t[d*j])) end do; Graph: (cf. Jerri [3, p. 92])Q[1] := pointplot([seq([t[n], X[n][1]], n = 0..4)], color=black, labels=[t, "x1"], labelfont=[times, bold, 12]): Q[2] := pointplot([seq([t[n], X[n][1]], n = 0..4)], color=blue, style = line, thickness=1, labels=[t, "x1"], labelfont=[times, bold, 12]): display(Q[1], Q[2]); Applying Volterra to (6.1) with step size h = 0.1:Volterra(1, .1, 63, 10^(-5), 5); N:=n-1; (total number of steps of size h)d:=3; (In the following table, every value is displayed if d = 1; every other value if d = 2, etc.) M:=21; (M + 1 is the number of values that will be displayed in each column.) Comparison of the numerical solution x1(t) with the analytical solution sin t: printf(" t x1(t) sin(t) \n"); for j from 0 to M do printf("%10.1f %10.4f %10.4f\n", t[d*j], X[d*j][1], sin(t[d*j])) end do; Graphs: x1(t) vs. t:Q[3] := pointplot([seq([t[n], X[n][1]], n = 0..63)], color=black, labels=[t, "x1"], labelfont=[times, bold, 12], symbol=box, symbolsize=6, legend="numerical solution"): display(Q[3]); sin t vs t:Q[4] := pointplot([seq([t[n], sin(t[n])], n = 0..63)], color=blue, labels=[t, "x"], labelfont=[times, bold, 12], symbol=cross, symbolsize=10, legend="x = sin(t)"): display(Q[4]); x1(t) and sin t plotted on the same set of axes: display(Q[3], Q[4], title="A comparision of the numerical and analytical solutions", titlefont=[times, bold, 12], labels=[t, " "]);

The unique analytical solution of the scalar Volterra integro-differential equation (7.1) NiMtSSVkaWZmR0kqcHJvdGVjdGVkR0YlNiQtSSJ4RzYiNiNJInRHRilGKw== = NiMsJComIiIkIiIiLUkieEc2IjYjSSJ0R0YpRiYhIiI= - NiMtSSRJbnRHNiI2JComLUkkZXhwR0YlNiMsJkkic0dGJSIiIkkidEdGJSEiIkYtLUkieEdGJTYjRixGLS9GLDsiIiFGLg==satisfying the initial condition x(0) = 1 is x(t) = e-2t - te-2t (cf. Burton [2, pp. 8-9]). In order to solve (7.1) with the Volterra procedure, we must first change (7.1) to an equivalent system of two Volterra integral equations by letting x1(t) = x(t) and x2(t) = x'(t). Then we obtain NiMvLSZJInhHNiI2IyIiIjYjSSJ0R0YnLCZGKUYpLUkkSW50RzYkSSpwcm90ZWN0ZWRHRjBJKF9zeXNsaWJHRic2JC0mRiY2IyIiIzYjSSJzR0YnL0Y4OyIiIUYrRik=(7.2) NiMtJkkieEc2IjYjIiIjNiNJInRHRiY= = NiMsJCIiJCEiIg== + NiMtSSRJbnRHNiRJKnByb3RlY3RlZEdGJkkoX3N5c2xpYkc2IjYkPCMsJiomLUkkZXhwR0YlNiMsJkkic0dGKCIiIkkidEdGKCEiIkYyLSZJInhHRig2I0YyNiNGMUYyRjQqJiIiJEYyLSZGNzYjIiIjRjlGMkY0L0YxOyIiIUYz. Since (7.2) consists of two equations, m = 2. m := 2; Enter the function f(t):f := proc(t) Vector([1, -3]) end proc: Check that f has been entered correctly:f(t);Enter the kernel k(t, s, x):k := proc(t, s, x::Vector) Vector([x[2], -exp(s-t)*x[1]-3*x[2]]); end proc:Check that k(t, s, x) has been entered correctly:v := Vector(m, symbol = x): vlist := convert(v, 'list'): k(t, s, v);Check that the system (7.2) has been entered correctly:for i from 1 to m do x[i](t) = f(t)[i] + Int(k(t, s, v)[i], s = 0..t) end do; Applying Volterra to (7.2) with step size h = 0.1:Volterra(2, .1, 40, 10^(-8), 5);N:=n-1; (total number of steps of size h)d:=2; (In the following table, every value is displayed if d = 1; every other value if d = 2, etc.) M:=20; (M + 1 is the number of values that will be displayed in each column.) Comparison of the exact values of the analytical solution x(t) = e-2t - te-2t with the numerical approximations of x(t): printf(" t exact approximation \n"); for j from 0 to M do printf("%10.1f %10.5f %10.5f\n", t[d*j], exp(-2*d*t[j])-d*t[j]*exp(-2*d*t[j]), X[d*j][1]) end do; Graphs: Graph of x1(t), the numerical approximation of x(t):R[1] := pointplot([seq([t[n], X[n][1]], n = 0..40)], color=black, labels=[t, "x1"], labelfont=[times, bold, 12]): R[2] := pointplot([seq([t[n], X[n][1]], n = 0..40)], color=blue, style = line, connect=true, thickness=1, labels=[t, "x1"], labelfont=[times, bold, 12], title="Numerical solution", titlefont=[times, bold, 12]): display(R[1], R[2]); Graph of the analytical solution x(t) = e-2t - te-2t:plot(exp(-2*t) - t*exp(-2*t), t=0..4, labels=[t, x], labelfont=[times, bold, 12], color=magenta, title="Analytical solution", titlefont=[times, bold, 12], thickness=2);

In this section, we numerically solve the Volterra system of three integral equations NiMtJkkieEc2IjYjIiIiNiNJInRHRiY= = NiMqJiIiKSIiIi1JJHNpbkc2JEkqcHJvdGVjdGVkR0YpSShfc3lzbGliRzYiNiNJInRHRitGJQ== + NiMtSSRJbnRHNiRJKnByb3RlY3RlZEdGJkkoX3N5c2xpYkc2IjYkKigtJkkieEdGKDYjIiIjNiNJInNHRigiIiItJkYtNiMiIiRGMEYyLUkkY29zR0YlRjBGMi9GMTsiIiFJInRHRig=(8.1) NiMtJkkieEc2IjYjIiIjNiNJInRHRiY= = NiMqJkkidEc2IiIiIiIiJSEiIg== + NiMtSSRpbnRHNiI2JDwjLCYqKi1JJGV4cEc2JEkqcHJvdGVjdGVkR0YtSShfc3lzbGliR0YlNiMsJkkic0dGJSIiIkkidEdGJSEiIkYyLUkkc2luR0YsNiNGM0YyLUkkY29zR0YsNiMqJiIiI0YyRjFGMkYyLSZJInhHRiU2I0YyNiNGMUYyRjIqJiIiJEYyLSZGPzYjRjxGQUYyRjQvRjE7IiIhRjM= NiMtJkkieEc2IjYjIiIkNiNJInRHRiY= = NiMsJiIiIkYkLUkkY29zRzYiNiNJInRHRichIiI= + NiMtSSRJbnRHNiRJKnByb3RlY3RlZEdGJkkoX3N5c2xpYkc2IjYkPCMsKComLSZJInhHRig2IyIiIjYjSSJzR0YoRjEtJkYvNiMiIiNGMkYxRjEqKi1JJGV4cEdGJTYjLCZGM0YxSSJ0R0YoISIiRjEtSSRzaW5HRiU2I0Y9RjEtSSRjb3NHRiVGMkYxRjRGMUYxLSZGLzYjIiIkRjJGPi9GMzsiIiFGPQ==and graph the results. Since (8.1) consists of three equations, m = 3.m:=3; Enter the function f(t):f := proc(t) Vector([8*sin(t), t/4, 1-cos(t)]) end proc: Check that f has been correctly entered:f(t); Enter the kernel k(t, s, x):k := proc(t, s, x::Vector) Vector([x[2]*x[3]*cos(s), exp(s-t)*sin(t)*cos(2*s)*x[1]-3*x[2], x[1]*x[2] + exp(s-t)*sin(t)*cos(s)*x[2] -x[3]]); end proc: Check that k(t, s, x) has been correctly entered:v := Vector(m, symbol = x): vlist := convert(v, 'list'): k(t, s, v);Check that (8.1) has been correctly entered:for i from 1 to m do x[i](t) = f(t)[i] + Int(k(t, s, v)[i], s = 0..t) end do; Applying Volterra to (8.1) with step size h = 0.1:Volterra(3, 0.1, 150, 10^(-8),5);N:=n-1; (total number of steps of size h)d:=5; (In the following table, every value is displayed if d = 1; every other value if d = 2, etc.) M:=30; (M + 1 is the number of values that will be displayed in each column.) Numerical approximation of the solution (x1(t), x2(t), x3(t)): printf(" t x1(t) x2(t) x3(t) \n"); for j from 0 to M do printf("%6.2f %8.4f %8.4f %8.4f\n", t[d*j], X[d*j][1], X[d*j][2], X[d*j][3]) end do; Graphs: x1(t) vs. t: S[1] := pointplot([seq([t[n], X[n][1]], n = 0..N)], color=black, labels=[t, "x1"], labelfont=[times, bold, 12]): S[2] := pointplot([seq([t[n], X[n][1]], n = 0..N)], color=green, style=line, thickness=1, labels=[t, "x1"], labelfont=[times, bold, 12]): display(S[1], S[2]); x2(t) vs. t: S[3] := pointplot([seq([t[n], X[n][2]], n = 0..N)], color=black, labels=[t, "x2"], labelfont=[times, bold, 12]): S[4] := pointplot([seq([t[n], X[n][2]], n = 0..N)], color=magenta, style=line, thickness=1, labels=[t, "x2"], labelfont=[times, bold, 12]): display(S[3], S[4]); x3(t) vs. t: S[5] := pointplot([seq([t[n], X[n][3]], n = 0..N)], color=black, labels=[t, "x3"], labelfont=[times, bold, 12]): S[6] := pointplot([seq([t[n], X[n][3]], n = 0..N)], color=gold, style=line, thickness=1, labels=[t, "x3"], labelfont=[times, bold, 12]): display(S[5], S[6]); Graph of the numerical approximation of the solution (x1(t), x2(t), x3(t)):S[7] := pointplot3d([seq([X[n][1], X[n][2], X[n][3]], n = 0..N)], color=black, labels=["x1", "x2", "x3"], labelfont=[times, bold, 12], tickmarks=[3,3,3], axes=normal, orientation=[50, 60]): S[8] := pointplot3d([seq([X[n][1], X[n][2], X[n][3]],n = 0..N)], color=blue, labels=["x1","x2","x3"],labelfont=[times, bold, 12], style=line, tickmarks = [3, 3, 3], thickness=1): display(S[7], S[8]);

This worksheet can be used to compute and graph numerical solutions of other systems of Volterra integral equations of the second kind by modifying the examples in Sections 5, 6, and 8. Numerical solutions of scalar Volterra integro-differential equations can also be computed and graphed by modifying the example in Section 7.

[1] Richard L. Burden and J. Douglas Faires, Numerical Analysis, 8th ed., Thomson Brooks/Cole, Belmont, CA, 2005.[2] T. A. Burton, Volterra Integral and Differential Equations, 2nd ed., Mathematics in Science & Engineering, 202, Elsevier, 2005.[3] Abdul J. Jerri, Introduction to Integral Equations with Applications, Pure & Applied Mathematics, vol. 93, Marcel Dekker, New York, 1985.[4] Peter Linz, Analytical and Numerical Methods for Volterra Equations, Studies in Applied Mathematics 7, SIAM, Philadelphia, 1985.[5] W. H. Press, S. A. Teukolsky, W. T. Vetterling, B. P. Flannery, Numerical Recipes in Fortran 77: the Art of Scientific Computing, vol. 1, 2nd ed., Cambridge University Press, 1986-1992. An Adobe Acrobat version is found at the web site: http://library.lanl.gov/numerical/ .

Legal Notice: The copyright for this application is owned by the authors. Neither Maplesoft nor the authors are responsible for any errors contained within and are not liable for any damages resulting from the use of this material. This application is intended for non-commercial, non-profit use only. Contact the authors for permission if you wish to use this application in for-profit activities.* End of worksheet *