MATH202: Differential Equations II

MAPLE Assignment Four

# linbvp.maple  Maple code to solve the linear boundary value problem
# 12.10.03      y_{xx} +q(x)*y = g(x)
#               y(x=0)    = ystart
#               y(x=xend) = yend 
#
with(linalg):
with(plots):

xstart := 0:
xend   := evalf(Pi):
ystart := 0:  # the solution at the gridpoint x=xstart
yend   := 0:  # the solution at the gridpoint x=xend
points := 20:
h      := (xend-xstart)/(points-1);

exactsol := x -> 32*(sin(x/2)+cos(x/2)-1) +4.0*(x**2-evalf(Pi**2)*sin(x/2));
p3 := plot(32*(sin(x/2)+cos(x/2)-1) +4*(x**2-Pi**2*sin(x/2)),x=xstart..xend,colour=RED):
q := x -> 0.25;
g := x -> x**2;

xpos := proc(n)            # calculate the position of the grid points.
        option remember;
        xpos(n-1)+h 
        end:
xpos(1) := 0:

xposition := [seq([xpos(n)],n=1..points)]:

# define the matrix A 
     # first put zeroes in every location
A := matrix(points,points,0):
     # define the non-zero numbers in the matrix

for i from 1 by 1 to points
do
   if i=1 then 
        A[i,i] := 1;
   elif i>1 and i<points then
        A[i,i-1] := 1;
        A[i,i]   := -2 +(h^2)*q(xpos(i));
        A[i,i+1] := 1;
   else
        A[i,i] := 1; 
   fi:
od:   

# define the vector c
c := matrix(points,1):
c[1,1] := ystart:
c[points,1] := yend:
for i from 2 by 1 to (points-1)
do
      c[i,1] := (h^2)*g(xpos(i));
od:


#print(xposition);
#print(A);
#print(c);

# Solve the matrix equation Ax=c
ysolution := linsolve(A,c);

# manipulate xposition and ysolution so we can plot the data

xdata := convert(transpose(convert(xposition,vector)),vector);
ydata := convert(ysolution,vector);

p1 := pointplot({seq([xdata[i],ydata[i]],i=1..points)},colour=BLACK):
p2 := listplot([seq([xdata[i],ydata[i]],i=1..points)],colour=BLACK):
display({p1,p2},title="linear bvp");
display({p1,p2,p3},title="Comparison of exact and numerical solution");

# calculate the relative percentage error at each interior point

seq(100.0*(exactsol(xdata[i])-ydata[i])/exactsol(xdata[i]),i=2..(points-1));

# calculate the absolute error at each interior point

AbsError := [seq(abs(exactsol(xdata[i])-ydata[i]),i=2..(points-1))];

# Find the maximum absolute error over all interior points

max(seq(abs(exactsol(xdata[i])-ydata[i]),i=2..(points-1)));


A      := 'A':
c      := 'c':
exactsol := 'exactsol':
g      := 'g':
h      := 'h':
p1     := 'p1':
p2     := 'p2':
p3     := 'p3':
points := 'points':
q      := 'q':
xdata  := 'xdata':
xend   := 'xend':
xpos   := 'xpos':
xstart := 'xstart':
ydata  := 'ydata':
yend   := 'yend':
ysolution := 'ysolution':
ystart := 'ystart':