MATH 542 NUMERICAL SOLUTIONS OF DIFFERENTIAL EQUATIONS

INTRODUCTION
 GRADIENT FIELDS WITH SOLUTIONS AUTONOMOUS SYSTEMS GRAD FIELD FOR A 2 EQ. SYSTEM WITH SOL FIRST ORDER PDE EXAMPLE TRAFFIC MODEL W/O SHOCK TRAFFIC MODEL WITH SHOCK
INITIAL VALUE PROBLEMS(ODE)
 INITIAL VALUE PROBLEMS (MAIN AND FUNCTION) RUNGE-KUTTA VARIABLE STEP RUNGE-KUTTA (1,2) BACKWARDS DIFFERENTIATION FORMULAS QUASI-VARIABLE STEP 2ND ORDER BDF ADAMS METHODS A-STABILITY REGION
FINITE DIFFERENCE METHODS
 EXAMPLE 1 EXAMPLE 2 HANGING CABLE EXAMPLE LAPLACE'S EQUATION HEAT EQUATION WAVE  EQUATION
FINITE ELEMENT METHOD
 LINEAR SUBDOMAIN COLLOCATION GALERKIN TRIANGULAR GALERKIN TRIANGULAR GRID REFINEMENT
HOME

Also requires the following function subprogram

function F=Ode_function(x,y)
global func_f
F=eval(func_f);

MAIN PROGRAM

%plots the direction field and solutions for the
%equation dy/dx = f(x,y)
%enter f as a character string MATLAB code for f
clear
global
func_f
clf reset

%enter grid limits and size from the keyboard
%x_min=input('x minimum = '); x_max=input('x maximum = ');
%y_min=input('y minimum = '); y_max=input('y maximum = ');
%delta_x=input('delta x = '); delta_y=input('delta y = ');
%or just edit the following line

x_min = 0; x_max = 3; y_min = 0; y_max = 4; delta_x = .2; delta_y = .2;
[x,y] = meshgrid(x_min:delta_x:x_max, y_min:delta_y:y_max);
[n,m]=size(x);

%enter righthand side function
%func_f=input('dy/dx = ','s');
%or just edit the following line

func_f ='-9*(x-1).*y'
Heading = ['\n\n\n', 'Window Parameters for the Equation \n\n',' dy/dx = ', func_f, '\n\n'];
fprintf('\n x-min = %g x-max = %g y-min = %g y-max = %g \n',...
x_min, x_max, y_min, y_max)
fprintf('\n delta-x = %g delta-y = %g \n\n', delta_x, delta_y)
z=eval(func_f);
One=ones(n,m);
L=sqrt(One.^2+z.^2);
quiver(x, y, One./L,z./L,.9)
axis([x_min,x_max,y_min,y_max])
grid on
xlabel('x axis')
ylabel('y axis')
Title=['dy/dx = ', func_f];
title(Title)
Test=input('Plot a solution? ', 's');

while Test(1) == 'y'

fprintf('\n Use the mouse to locate the initial values in the Figure Window \n')
[int_x, int_y]=ginput(1);
hold on
plot(int_x,int_y,'ko','LineWidth',2,'MarkerSize',7)
[u,v]=ode45('Ode_function',[int_x, x_max], int_y);
plot(u,v, 'r','LineWidth',2)
Test = input('Plot another solution? ', 's');

end

 TOP BACK

AUTONOMOUS SYSTEMS
Also requires the following function subprogram

function F=ASystem_function(t,X)
global func_f func_g
x=X(1);
y=X(2);
F=[eval(func_f); eval(func_g)];

MAIN PROGRAM

%plots the direction field for the
%autonomous system dx/dt = f(x,y) dy/dt= g(x,y)
%enter f and g as a character string MATLAB code for f

clear
global func_f func_g
clf reset

%enter grid limits and size from the keyboard
%x_min=input('x minimum = '); x_max=input('x maximum = ');
%y_min=input('y minimum = '); y_max=input('y maximum = ');
%delta_x=input('delta x = '); delta_y=input('delta y = ');
%or just edit the following line

x_min = -3; x_max = 3; y_min = -3; y_max = 3; t_max=10; delta_x = .2; delta_y = .2;
[x,y] = meshgrid(x_min:delta_x:x_max, y_min:delta_y:y_max);
[n,m]=size(x);

%enter righthand side function
%func_f=input('dx/dt = ','s');
%func_g=input('dy/dt = ','s');
%or just edit the following lines

func_g ='-sin(x)';
func_f ='y';
Heading = ['\n\n\n', 'Window Parameters for the Autonomous System \n\n',' dx/dt = ',...
func_f, ' dy/dt = ', func_g,'\n\n'];
fprintf('\n x-min = %g x-max = %g y-min = %g y-max = %g \n',...
x_min, x_max, y_min, y_max)
fprintf('\n delta-x = %g delta-y = %g \n\n', delta_x, delta_y)
fprintf('\n t-max = %g \n\n', t_max)
zn=eval(func_g);
zd=eval(func_f);
L=sqrt(zn.^2+zd.^2);
quiver(x, y, zd./L,zn./L,.7)
axis([x_min,x_max,y_min,y_max])
grid on
xlabel('x axis')
ylabel('y axis')
Title=['dx/dt = ', func_f,' dy/dt = ', func_g];
title(Title)
Test=input('Plot a solution? ', 's');

while Test(1) == 'y'

fprintf('\n Use the mouse to locate the initial values in the Figure Window \n')
[int_x, int_y]=ginput(1);

for k=1:t_max
[t,V]=ode45('ASystem_function',[0,1], [int_x,int_y]);
hold on
plot(int_x,int_y,'ko','LineWidth',2)
NumberString=[' t = ',int2str(k-1)];
text(int_x,int_y, NumberString,'fontsize',12,'fontweight','bold','horizontalalignment','left')
plot(V(:,1),V(:,2), 'r','LineWidth',2)
nn=length(V(:,1));
int_x=V(nn,1); int_y=V(nn,2);

end

plot(int_x,int_y,'ko','LineWidth',2)
NumberString=[' t = ',int2str(k)];
text(int_x,int_y, NumberString,'fontsize',12,'fontweight','bold','horizontalalignment','left')
Test = input('Plot another solution? ', 's');

end

 TOP BACK

GRADIENT FIELD FOR A SYSTEM OF 2 EQUATIONS WITH SOLUTION

Also requires this function subprogram

function F=OdeSys_function(x,u)

global func_f func_g
y=u(1); z=u(2);
F=[eval(func_f)
eval(func_g)];

Main Program

%plots the direction field and solutions for the

%system dy/dx = f(x,y,z) dz/dx = g(x,y,z)

%enter f and g as a character string

clear

global func_f func_g

clf reset

%set the figure window size and grid size

x_min = 0; x_max = 20; y_min = -2; y_max = 2; z_min = -2; z_max = 2;

delta_x = 2; delta_y = 1; delta_z = 1;

[x,y,z] = meshgrid(x_min:delta_x:x_max, y_min:delta_y:y_max, z_min:delta_z:z_max);

[n,m,p]=size(x);

%enter righthand side functions f and g

func_f='z'

func_g ='-sin(y)';

%enter the initial conditions at x=0

int_y=0; int_z=1.5;

Heading = ['\n\n\n', 'Window Parameters for the System \n\n',' dy/dx = ', func_f, '\n\n',' dz/dx = ', func_g, '\n\n'];

fprintf('\n x-min = %g x-max = %g y-min = %g y-max = %g z-min = %g z-max = %g \n',...

x_min, x_max, y_min, y_max, z_min, z_max)

fprintf('\n delta-x = %g delta-y = %g delta-z = %g \n\n', delta_x, delta_y, delta_z)

u=eval(func_f);

v=eval(func_g);

One=ones(n,m,p);

L=sqrt(One.^2+u.^2+v.^2);

quiver3(x, y, z, One./L,u./L,v./L,.9);

axis([x_min,x_max,y_min,y_max,z_min,z_max])

grid on

%axis labels

xlabel('time axis')

ylabel('angle axis')

zlabel('velocity axis')

Title=['dy/dx = ', func_f,' dz/dx = ', func_g];

title(Title)

rotate3d on, hold on

plot3(0,int_y,int_z,'ok','LineWidth',2, 'MarkerSize',7)

[uu,vv]=ode45('OdeSys_function',[0, x_max], [int_y; int_z]);

plot3(uu,vv(:,1),vv(:,2), 'r','LineWidth',2)

 TOP BACK

FIRST ODER PDE EXAMPLE

%solution of a first order PDE

clear, clf reset
[x,t]=meshgrid(0:.2:24,0:.1:2);
u=cos(x).*exp(-x.*t.*t/2);
rotate3d on
surf(x,t,u,t)
title('Graph of u = cos(x) exp(-x t^2/2)')
xlabel('x axis'),ylabel('t axis'),zlabel('u axis')

 TOP BACK

TRAFFIC MODEL W/O SHOCK

%Traffic Model without shock

clear
clf
x = -3:.1:3;
t = 0:.25:4;

%calculation of d0 = d(x,0)

xx = -1:.1:1;
w = 1/3 + 1/12*(xx+1).*(xx+1).*(2-xx);
d0 = [1/3*ones(1,20),w,2/3*ones(1,20)];
n=length(x);
m=length(t);

for i=1:n,

for j=1:m,

X(i,j)=(1-2*d0(i))*t(j)+x(i);
T(i,j)=t(j);
D(i,j)=d0(i);

end

end

surf(X,T,D,X)
colormap hsv
view(-20,30)
colorbar
grid on
xlabel('x axis')
ylabel('t axis')
zlabel('d axis')
rotate3d on

 TOP BACK

TRAFFIC MODEL WITH SHOCK

%Traffic Model with shock

clear
x = -3:.1:3;
t = 0:.25:4;

%calculation of d0 = d(x,0)

xx = -1:.1:1;
w = 1/3 + 1/12*(xx+1).*(xx+1).*(2-xx);
d0 = [1/3*ones(1,20),w,2/3*ones(1,20)];
n=length(x);
m=length(t);
for i=1:n,
for j=1:m,
X(i,j)=(1-2*d0(i))*t(j)+x(i);
T(i,j)=t(j);
D(i,j)=d0(i);
end
end

c=ones(n,m);
for j=1:m
a=0; b=0;
for i=1:n
if X(i,j)>0 & a==0
is=i; a=1;
end
if
X(n-i+1,j)<0 & b==0
ib=n-i+1; b=1;
end
end
if
is<ib
for k=is:ib
X(k,j)=0;
c(k,j)=0;
end
end
end

surf(X,T,D,X)
colormap hsv
view(-20,30)
colorbar
grid on
xlabel('x axis')
ylabel('t axis')
zlabel('d axis')
rotate3d on

 TOP BACK

EXAMPLE 1

%Example 1
clear, clf
C = 1; D = .01; K = .1; A=[]; n = 4;

for i = 1:8
flops(0)
h = 1/n;
%Construct the Tridiagonal System
l=ones(1,n-2); d=-(2+K*h^2/D)*ones(1,n-1); b=[zeros(1,n-2), -C];
%Solve the Tridiagonal System
U = tridiag(l,d,l,b);
p=flops;
%Plot the Approx. Solution
X = 0:h:1; U = [0, U, C];
plot(X,U,'b')
grid on, xlabel('x-axis'), ylabel('u-axis'), title('Example 2.6')
hold on
%Save values at the midpoint for error analysis
A=[A;n,p,X(n/2+1),U(n/2+1)];
n = n*2;
end

hold off
%error analysis
[nn,mm]=size(A);
diff(1)=0;ratio(1)=0; ratio(2)=0;

for i=1:nn-1,
diff(i+1)=A(i+1,4)-A(i,4);
end

for i=2:nn-1,
ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,4) + diff(i)/(1-ratio(i+1));
end

z=[A(:,1)';A(:,2)';A(:,4)'; diff; ratio; xacc];
fprintf(1,'\n\n n flops y diff ratio y acc\n\n')
fprintf(1,'%3.0f %7.0f %16.14f %13.10f %13.10f %16.14f\n',z)

 TOP BACK

EXAMPLE 2

%Example 2
clear, clf
C = 1; D = .01; K = .1; A=[];n = 4;
for i = 1:8
h = 1/n;
flops(0)
%Using Backward Diff at endpoint
%l=[ones(1,n-2),-1]; u=ones(1,n-1); d=[-(2+K*h^2/D)*ones(1,n-1),1];
%b=[zeros(1,n-1), h*C/D];Y = tridiag(l,d,u,b);
%Using an imaginary node at x=1+h

l=[ones(1,n-2),2]; u=ones(1,n-1); d=[-(2+K*h^2/D)*ones(1,n)];
b=[zeros(1,n-1), -2*h*C/D];Y = tridiag(l,d,u,b);
%Using a second order approx at x=1
%l=[ones(1,n-2),-4];u=ones(1,n-1);d=[-(2+K*h^2/D)*ones(1,n-1),3];
%b=[zeros(1,n-1), 2*h*C/D];
%M=zeros(n,n)+diag(d)+diag(l,-1)+diag(u,1);M(n,n-2)=1;
%Y=(M\b')';

p=flops;
X = 0:h:1; Y = [0, Y];
plot(X,Y,'b')
grid on, xlabel('x-axis'), ylabel('u-axis'), title('Example 2.7')
hold on
%Save values at the midpoint for error analysis
A=[A;n,p,X(n/2+1),Y(n/2+1)];
n = n*2;
end
hold off
A;
[nn,mm]=size(A);
diff(1)=0;ratio(1)=0; ratio(2)=0;
for i=1:nn-1,
diff(i+1)=A(i+1,4)-A(i,4);
end
for
i=2:nn-1,
ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,4) + diff(i)/(1-ratio(i+1));
end
z=[A(:,1)';A(:,2)';A(:,4)'; diff; ratio; xacc];
fprintf(1,'\n\n n flops u(0.5) diff ratio u acc\n\n')
fprintf(1,'%3.0f %7.0f %16.14f %13.10f %13.10f %16.14f\n',z)

 TOP BACK

HANGING CABLE EXAMPLE

This code calculates the shape of a cable hanging under its own weight and then provides a graph.  The non-linear equations are solved by Newton's Method.

MAIN PROGRAM

clear, clf
L=3; a=1; b=2; Ha=1; Hb=2;s=-5;A=[];n=4;
for ii=1:8
flops(0)
%Initial Guess for Newton's Method
h=(b-a)/n; x=a-h:h:b+h;
y=(Hb-Ha)/((b-a)^2)*(x-a).^2+Ha+(s-(s/(b-a))*(x-a)).*(x-a);
%Newton's Method to solve the system
for i=1:8
r=DFHC(y,Ha,Hb,h,L)\(-FHC(y,Ha,Hb,h,L));
y=y+r';
end
p=flops;
xx=x(2:n+2);
yy=y(2:n+2);
plot(xx,yy)
grid on,hold on, xlabel('x axis'), ylabel('y axis')
A=[A;n,p,x(n/2+2),y(n/2+2)];
n=2*n;
end
[nn,mm]=size(A);
diff(1)=0;ratio(1)=0; ratio(2)=0;
for i=1:nn-1,
diff(i+1)=A(i+1,4)-A(i,4);
end
for
i=2:nn-1,
ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,4) + diff(i)/(1-ratio(i+1));
end
z=[A(:,1)';A(:,2)';A(:,4)'; diff; ratio; xacc];
fprintf(1,'\n\n n flops y(1.5) diff ratio y acc\n\n')
fprintf(1,'%3.0f %8.0f %16.14f %13.10f %13.10f %16.14f\n',z)

FUNCTION SUBPROGRAM

function z=FHC(y,Ha,Hb,h,L)
m=length(y);
Lc=1/(4*L*L); hc=1/(h*h);

for i=1:m-2
z(i)=hc*(y(i+2)-2*y(i+1)+y(i))^2-Lc*((y(m)-y(m-2)-y(3)+y(1))^2)*(1+(hc/4)*(y(i+2)-y(i))^2);
end

z(m-1)=y(2)-Ha;
z(m)=y(m-1)-Hb;
z=z';

JACOBIAN FUNCTION SUBPROGRAM

function M=DFHC(y,Ha,Hb,h,L)
m=length(y);
Lc=1/(4*L*L); hc=1/(h*h);

for i=1:m-2
d(i)=2*hc*(y(i+2)-2*y(i+1)+y(i))+Lc*((y(m)-y(m-2)-y(3)+y(1))^2)*(hc/2)*(y(i+2)-y(i));
u1(i)=-4*hc*(y(i+2)-2*y(i+1)+y(i));
u2(i)=2*hc*(y(i+2)-2*y(i+1)+y(i))-Lc*((y(m)-y(m-2)-y(3)+y(1))^2)*(hc/2)*(y(i+2)-y(i));
k1(i)=-2*Lc*(y(m)-y(m-2)-y(3)+y(1))*(1+(hc/4)*(y(i+2)-y(i))^2);
k2(i)=2*Lc*(y(m)-y(m-2)-y(3)+y(1))*(1+(hc/4)*(y(i+2)-y(i))^2);
k3(i)=2*Lc*(y(m)-y(m-2)-y(3)+y(1))*(1+(hc/4)*(y(i+2)-y(i))^2);
k4(i)=-2*Lc*(y(m)-y(m-2)-y(3)+y(1))*(1+(hc/4)*(y(i+2)-y(i))^2);
end

M=diag([d,0,0])+diag([u1,0],1)+diag(u2,2);

for i=1:m-2
M(i,1)=k1(i)+M(i,1); M(i,3)=k2(i)+M(i,3); M(i,m-2)=k3(i)+M(i,m-2); M(i,m)=k4(i)+M(i,m);
end

M(m-1,2)=1; M(m,m-1)=1;

APPROXIMATE JACOBIAN FUNCTION SUBPROGRAM

Use either the Jacobian Function or the Approximate Jacobian Function NOT BOTH

function J=DFHC(y,Ha,Hb,h,L)

%Computes an Approximate Jacobian for the given function F
%and evaluates at the vector y
%Assumes that the function 'f' will return a column vector

epsi=10^-7;
n=length(y);J=[];

for i=1:n

e=zeros(1,n); e(i)=epsi;
x=y+e; fip=FHC(x,Ha,Hb,h,L);
x=y-e; fim=FHC(x,Ha,Hb,h,L);
J=[J,(fip-fim)/(2*epsi)];

end

 TOP BACK

INITIAL VALUE PROBLEMS
Requires the following Main Program and two function subprograms.

Main
clear all
global func_f func_g
clf reset

%Initial step size h, starting and ending values of a<x<b
h=.125; a=0; b=5.3;
%enter righthand side function
func_f ='-9*(x-1).*y';
%func_f ='[y(2),sqrt(1+y(2)^2)]';
equation = ['dy/dx = ', func_f]
%enter the Jacobian of the righthand side
func_g = '-9*(x-1)'; ['Jacobian = ', func_g]
%func_g = '[0,1;0,y(2)/sqrt(1+y(2)^2)]'; ['Jacobian = ', func_g]
%Initial condition

ya=[1];
for j=1:6 flops(0)
[x,y] = feuler(a,b,h,ya);
%choose the method
p=flops;
n=length(x);
A(j,3)=y(n,1); A(j,2)=p; A(j,1)=h;
h=h/2;
plot(x,y(:,1),
'b')
grid on
xlabel(
'x-axis')
ylabel(
'y-axis')
title(equation)
hold on

end
%error estimation
hold off
[nn,mm]=size(A);diff=zeros(1,nn); ratio=zeros(1,nn); xacc=zeros(1,nn);

for i=1:nn-1,
diff(i+1)=A(i+1,3)-A(i,3);

end
for
i=2:nn-1,
ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,3) + diff(i)/(1-ratio(i+1));

end
z=[A(:,1)';A(:,2)';A(:,3)'; diff; ratio; xacc];
fprintf(1,
'\n\n h flops y diff ratio y acc\n\n')
fprintf(1,
'%7.5f %7.0f %16.14f %13.10f %13.10f %16.14f\n',z)

Required Function Subprograms

function F=Ode_function(x,y)
global func_f
F=eval(func_f);

function F=Jacobian_function(x,y)
global func_g
F=eval(func_g);

 TOP BACK

RUNGE-KUTTA

Forward Euler

function [u,v] = feuler(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
for i=1:n-1,
v(i+1,:) = v(i,:) + h*Ode_function(u(i),v(i,:));
end

Heun

function [u,v] = heun(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
for i=1:n-1,
w= v(i,:) + h*Ode_function(u(i),v(i,:));
v(i+1,:) = v(i,:)+ (.5)*h*(Ode_function(u(i),v(i,:))+Ode_function(u(i+1),w));
end

4th Order RK

function [u,v] = rk4(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
for i=1:n-1,
u0=u(i); h2=h/2;
v0=v(i,:); f0=Ode_function(u0,v0);
v1=v0+h2*f0; f1=Ode_function(u0+h2,v1);
v2=v0+h2*f1; f2=Ode_function(u0+h2,v2);
v3=v0+h*f2; f3=Ode_function(u0+h,v3);
v(i+1,:) = v0+(h/6)*(f0+2*f1+2*f2+f3);
end

4th Order RKF

```function [u,v] = rkf4(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
for i=1:n-1,
u0=u(i);
v0=v(i,:); f0=h*Ode_function(u0,v0);
v1=v0+f0/4; f1=h*Ode_function(u0+h/4,v1);
v2=v0+3*f0/32+9*f1/32; f2=h*Ode_function(u0+3*h/8,v2);
v3=v0+1932*f0/2197-7200*f1/2197+7296*f2/2197; f3=h*Ode_function(u0+12*h/13,v3);
v4=v0+439*f0/216-8*f1+3680*f2/513-845*f3/4104; f4=h*Ode_function(u0+h,v4);
v(i+1,:) = v0+(25*f0/216+1408*f2/2565+2197*f3/4104-f4/5);
end```

 TOP BACK

BACKWARDS DIFFERENTIATION FORMULA

All of these formulas require the Jacobian.

1st Order Gear or BDF Method

function [u,v] = gear1(a,b,h,ya)
u= (a:h:b)';
n=length(u); m=length(ya);
v=ya(:)';
for i=1:n-1,
p=v(i,:)+h*Ode_function(u(i),v(i,:));
for j =1:2
c = v(i,:)+h*Ode_function(u(i+1),p);
J = h*Jacobian_function(u(i+1),p) - eye(m);
V = J\(p-c)';
v(i+1,:) = p + V';
p=v(i+1,:);
end
end

2nd Order BDF Method

```function [u,v] = BDF2(a,b,h,ya)
u= (a:h:b)';
n=length(u); m=length(ya);
%Runge-Kutta Start
[uu,v] = heun(a,a+2*h,h,ya);

%Compute Beginning Diffs
D1(1,:)=(v(2,:)-v(1,:));D1(2,:)=(v(3,:)-v(2,:));
D2(1,:)=(D1(2,:)-D1(1,:));
for i=3:n-1,
p= v(i,:)+D1(i-1,:)+D2(i-2,:);
t=1; j=1;
while (t > 1e-16) & (j<6)
c= (4/3)*v(i,:)-(1/3)*v(i-1,:)+(2/3)*h*Ode_function(u(i+1),p);
J=2/3*h*Jacobian_function(u(i+1),p) - eye(m);
p1= (J\(p-c)')' + p;
t= norm(p1-p,'inf');
j=j+1; p=p1;
end
v(i+1,:) = p;
%Update the Diffs
D1(i,:)=(p-v(i,:)); D2(i-1,:)=(D1(i,:)-D1(i-1,:));
end```

4th Order BDF Method

```function [u,v] = BDF4(a,b,h,ya)
u= (a:h:b)';
n=length(u); m=length(ya);
%Runge-Kutta Start
[uu,v]=rk4(a,a+4*h,h,ya);

%Compute the Beginning Diffs
D1(1,:)=(v(2,:)-v(1,:));D1(2,:)=(v(3,:)-v(2,:));
D1(3,:)=(v(4,:)-v(3,:));D1(4,:)=(v(5,:)-v(4,:));
D2(1,:)=(D1(2,:)-D1(1,:));D2(2,:)=(D1(3,:)-D1(2,:));
D2(3,:)=(D1(4,:)-D1(3,:));
D3(1,:)=(D2(2,:)-D2(1,:));D3(2,:)=(D2(3,:)-D2(2,:));
D4(1,:)=(D3(2,:)-D3(1,:));

for i=5:n-1,
p=v(i,:)+D1(i-1,:)+D2(i-2,:)+D3(i-3,:)+D4(i-4,:);
t=1; j=1;
while (t > 1e-16) & (j<6)
c=(1/25)*(48*v(i,:)-36*v(i-1,:)+16*v(i-2,:)-3*v(i-3,:)+12*h*Ode_function(u(i+1),p));
J = 12/25*h*Jacobian_function(u(i+1),p) - eye(m);
p1= (J\(p-c)')' + p;
t= norm(p1-p,'inf');
j=j+1; p=p1;
end
v(i+1,:) = p;
%Update the Diffs
D1(i,:)=(p-v(i,:)); D2(i-1,:)=(D1(i,:)-D1(i-1,:));
D3(i-2,:)=(D2(i-1,:)-D2(i-2,:));D4(i-3,:)=(D3(i-2,:)-D3(i-3,:));
end```

 TOP BACK

VARIABLE STEP RUNGE-KUTTA (1, 2)

```%Runge Kutta 1-2 with automatic step size selection and graph of step size
clear all, clf reset
global func_f
%starting and ending values of a<x<b and local error tolerance
a=0; b=25; tol=.01;
%enter righthand side function
%func_f ='-9*(x-1).*y';
func_f='[y(2),-sin(y(1))]';
equation = ['dy/dx = ', func_f]
%Initial condition
ya=[1,0];
xout(1)=a; yout(1,:)=ya; hmax=(b-a)/20;
hmin=hmax*10^(-4);
x=xout(1); y=yout(1,:);
%choose a starting step size using the Shampine rule
p=norm(eval(func_f),'inf');
if p>0
sout(1)=min(hmax,sqrt(tol*max(norm(y,'inf'),1.0))/p);
else
sout(1)=hmax/10;
end
h=sout(1)
count=1;
while (x<b) & (h>hmin)& (count<1600)
count=count+1;
%set last point to the end point
if x+h>b
h=b-x;
end
%calculate the RK steps
k1=Ode_function(x,y);
k2=Ode_function(x+h,y+h*k1);
%calculate the error estimate
errest=norm(h/2*(k1-k2),'inf');
%use rel err if |y|>1 and abs err if |y|<1
tolerr=tol*max(norm(y,'inf'),1.0);
if errest<=tolerr
x=x+h; xout=[xout,x];
y=y+h/2*(k1+k2); yout=[yout;y];
sout=[sout,h];
end
%calculate new step
if errest~=0.0
h=min(hmax,0.9*h*sqrt(tolerr/errest));
end
end
if (x<b)
display('Integration did not reach end point requested')
x
end
[count, h, x, y, errest, tolerr]
subplot(2,1,1)
plot(xout,yout(:,1)')
title(equation)
grid on, xlabel('x-axis'), ylabel('y-axis')
subplot(2,1,2)
plot(xout,sout)
grid on, xlabel('x-axis'), ylabel('step size')```

 TOP BACK

QUASI-VARIABLE STEP 2ND ORDER BDF

```%BDF2 with automatic step size selection and graph of step size
clear all, clf reset
global func_f  func_g
%starting and ending values of a<x<b and local error tolerance
a=0; b=25; tol=.01;
%enter righthand side function and initial condition
%func_f ='-9*(x-1).*y'; ya=[1];
func_f='[y(2),-sin(y(1))]'; ya=[1,0];
equation = ['dy/dx = ', func_f]
%enter the Jacobian of the righthand side
%func_g = '-9*(x-1)';
func_g = '[0, 1;-cos(y(1)),0]';
['Jacobian = ', func_g]
m=length(ya);
xout(1)=a; yout(1,:)=ya; hmax=(b-a)/2;
hmin=hmax*10^(-10);
x=xout(1); y=yout(1,:);
%choose a starting step size
h=hmax/100; sout(1)=h;```
```%Use Second Order RK to obtain starting values
%calculate the RK steps
k1=Ode_function(x,y);
k2=Ode_function(x+h,y+h*k1);
%calculate the error estimate
errest=norm(h/2*(k1-k2),'inf');
%use rel err if |y|>1 and abs err if |y|<1
tolerr=tol*max(norm(y,'inf'),1.0);
%calculate starting step size
if errest~=0.0
h=min(hmax,0.9*h*sqrt(tolerr/errest));
end
h=h/2;
%take two steps
for i=1:2
k1=Ode_function(x,y);
k2=Ode_function(x+h,y+h*k1);
x=x+h; xout=[xout,x];
y=y+h/2*(k1+k2); yout=[yout;y];
sout=[sout,h];
end
%Compute Beginning Diffs
X=xout; Y=yout;
D1(1,:)=(Y(2,:)-Y(1,:));D1(2,:)=(Y(3,:)-Y(2,:));
D2(1,:)=(D1(2,:)-D1(1,:));
count=1; index=3;```
```%main loop
while (X(index)<b) & (h>hmin)& (count<7500)
count=count+1;
%calculate the BDF2 step
pp= Y(index,:)+D1(index-1,:)+D2(index-2,:);
t=1; j=1; p=pp;
while (t > 1e-16) & (j<6)
c= (4/3)*Y(index,:)-(1/3)*Y(index-1,:)+(2/3)*h*Ode_function(X(index)+h,p);
J=2/3*h*Jacobian_function(X(index)+h,p) - eye(m);
p1= (J\(p-c)')' + p;
t= norm(p1-p,'inf');
j=j+1; p=p1;
end```
```    %calculate the error estimate
errest=(2/11)*norm((pp-p),'inf');
```
```    %use rel err if |y|>1 and abs err if |y|<1
tolerr=tol*max(norm(p,'inf'),1.0);
%take step if true
if errest<=tolerr
X(index+1)=X(index)+h; xout=[xout,X(index+1)];
Y(index+1,:)=p; yout=[yout;Y(index+1,:)]; sout=[sout,h];
%calculate possible new step
hnew=(tolerr/errest)^(1/3)*h;
if hnew>2*h & index>4 & 2*h < hmax %remesh at 2h if true
h=2*h;
X(1)=X(index-4);X(2)=X(index-2);X(3)=X(index);
Y(1,:)=Y(index-4,:);Y(2,:)=Y(index-2,:);Y(3,:)=Y(index,:);
D1(1,:)=(Y(2,:)-Y(1,:));D1(2,:)=(Y(3,:)-Y(2,:));
D2(1,:)=(D1(2,:)-D1(1,:)); index=3;
else
D1(index,:)=(Y(index+1,:)-Y(index-1,:)); D2(index-1,:)=(D1(index,:)-D1(index-1,:));
index=index+1;
end
else
%remesh at h/2
h=h/2;
X(1)=X(index-1);X(2)=X(1)+h;X(3)=X(index);
Y(1,:)=Y(index-1,:);Y(2,:)=Y(index,:)-D1(index-1,:)/2-D2(index-2,:)/8;Y(3,:)=Y(index,:);
D1(1,:)=(Y(2,:)-Y(1,:));D1(2,:)=(Y(3,:)-Y(2,:));
D2(1,:)=(D1(2,:)-D1(1,:)); index=3;
end
end %end of main loop```
```if (X(index)<b)
display('Integration did not reach end point requested')
X(index)
end```
```[count, h, index, errest, tolerr]
subplot(2,1,1)
plot(xout,yout(:,1)')
%plot(yout(:,1),yout(:,2))
title(equation)
grid on, xlabel('x-axis'), ylabel('y-axis')
subplot(2,1,2)
semilogy(xout,sout)
grid on, xlabel('x-axis'), ylabel('step size')```

 TOP BACK

function [u,v] = center(a,b,h,ya)
u= (a:h:b)';
n=length(u);
v=ya(:)';
v(2,:)=ya+h*Ode_function(a,ya);
for i=1:n-2,
v(i+2,:) = v(i,:)+ 2*h*Ode_function(u(i+1),v(i+1,:));
end

u= (a:h:b)';
n=length(u);
[uu,v] = rk4(a,a+3*h,h,ya);
dv(1,:)=Ode_function(u(1),v(1,:));
dv(2,:)=Ode_function(u(2),v(2,:));
dv(3,:)=Ode_function(u(3),v(3,:));
dv(4,:)=Ode_function(u(4),v(4,:));
for i=4:n-1,
v(i+1,:)=v(i,:)+(h/24)*(55*dv(i,:)-59*dv(i-1,:)+37*dv(i-2,:)-9*dv(i-3,:));
dv(i+1,:)=Ode_function(u(i+1),v(i+1,:));
v(i+1,:)=v(i,:)+(h/24)*(9*dv(i+1,:)+19*dv(i,:)-5*dv(i-1,:)+dv(i-2,:));
dv(i+1,:)=Ode_function(u(i+1),v(i+1,:));
end

 TOP BACK

A-STABILITY REGION CALCULATION

The first program is for Runge-Kutta and the second for BDF methods

%Computation of the A-Stability Region
clear, clf reset
x=linspace(-3,1.5,200);
y=linspace(-3.5,3.5,200);
[X,Y]=meshgrid(x,y);
Z=X +Y*i;
%Euler's Method
M=abs(1+Z);
[c,h]=contour(X,Y,M,[1,1]);
set(h,'linewidth',2,'edgecolor','g')
%colormap([0,0,1])
hold on
%Heun's Method
M=abs(1+Z/2.*(2+Z));
[c,h]=contour(X,Y,M,[1,1]);
set(h,'linewidth',2,'edgecolor','r')
%RK4
M=abs(1+Z/6.*(1+Z.*(1+Z/2.*(1+Z/2))));
[c,h]=contour(X,Y,M,[1,1]);
set(h,'linewidth',2,'edgecolor','k')
grid on
axis equal
title('Runge-Kutta A-Stabiity Regions')
%RKF4
M=abs(1+Z.*(1+Z/2.*(1+Z/3.*(1+Z/4.*(1+Z/(13/3))))));
[c,h]=contour(X,Y,M,[1,1]);
set(h,'linewidth',2,'edgecolor','b')
legend('Eulers Method', 'Heuns Method','RK4','RKF4')

```%Computation of the A-Stability Region(BDF1 and BDF2)
clear, clf reset
x=linspace(-3,3,21);
y=linspace(-3,3,21);
[X,Y]=meshgrid(x,y);
Z=X +Y*i;
%Backward Euler's Method
M=abs(1./(1-Z));
mesh(X,Y,M)
rotate3d on
colormap('copper')
axis([-3,3,-3,3,0,10])
title('Stability Surface for Backward Euler')
%BDF2 Method
M=abs((2+sqrt(1+2*Z))./(3-2*Z));
figure(2)
clf
mesh(X,Y,M)
rotate3d on
colormap('copper')
axis([-3,3,-3,3,0,10])
title('Stability Surface for BDF2 (first root)')
M=abs((2-sqrt(1+2*Z))./(3-2*Z));
figure(3)
clf
mesh(X,Y,M)
rotate3d on
colormap('copper')
axis([-3,3,-3,3,0,1])
title('Stability Surface for BDF2 (second root)')```

 TOP BACK

LAPLACE'S EQUATION

Class Example:  System Matrix is explicitly calculated

clear all
%Irregular Geometry Finite Diff Elliptic PDE Example
h=0.5; n=9; maxer=1; count=0;
%Molecule vertex list
M=[15,2,5,2;16,3,6,1;17,4,7,2;3,8,8,3;1,6,10,6;2,7,11,5;3,8,12,6;4,9,13,7;8,14,14,8]
A=zeros(17,17);B=zeros(17,1);

for i=1:9
A(i,i)=-4;

for j=1:4
A(i,M(i,j))=A(i,M(i,j))+1;

end
end

%Boundary Conditions
for i=10:13
A(i,i)=1;A(i,i-9)=-1;A(i,i-5)=2*h*0.2; B(i)=2*h*0.2*20;

end
A(14,14)=1;A(14,8)=-1;A(14,9)=2*h*0.2; B(14)=2*h*0.2*20;
for i=15:17
A(i,i)=1; B(i)=100;

end
fprintf('\n %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f \n',A')
B
u=(A\B)'
%plot the result
x=0:.5:2;
y=0:.5:1;
z=[100,100,100,u(3),u(7);u(1:4),u(8);u(5:9)]
hold off
surf(x,y,z)
view(160,20)
rotate3d on
colorbar

Code Using Gauss-Seidel

clear

%Irregular Geometry Finite Diff Elliptic PDE Example
h=0.5; n=9; maxer=1; count=0;

%Molecule vertex list
M=[15,2,5,2;16,3,6,1;17,4,7,2;3,8,8,3;1,6,10,6;2,7,11,5;3,8,12,6;4,9,13,7;8,14,14,8]

%Initialize solution
u=80*ones(1,17);
u(15)=100; u(16)=100; u(17)=100;

%Solve the linear system with Gauss-Seidel
flops(0);
while maxer > .00001
count=count+1;
maxer=0;

for i=1:9
r=(u(M(i,1))+u(M(i,2))+u(M(i,3))+u(M(i,4)))/4;
maxer=max(maxer,abs(r-u(i)));
u(i)=r;
end
for i=1:4
u(9+i)=2*h*(0.2)*(20-u(i+4))+u(i);
end
u(14)=2*h*(0.2)*(20-u(9))+u(8);
end
count, flops
u(1:9);

%plot the result
x=0:.5:2;
y=0:.5:1;
z=[100,100,100,u(3),u(7);u(1:4),u(8);u(5:9)]
hold off
surf(x,y,z)
view(160,20)
rotate3d on
colorbar

Using SOR

clear

%Irregular Geometry Finite Diff Elliptic PDE Example

h=0.25; n=30; maxer=1; count=0;w=1.65

%Molecule vertex list

M=[40 2 7 2;41 3 8 1;42 4 9 2;43 5 10 3; 44 6 11 4; 5 12 12 5;

1 8 14 8;2 9 15 7;3 10 16 8;4 11 17 9; 5 12 18 10;6 13 19 11;12 20 20 12;

7 15 22 15;8 16 23 14;9 17 24 15;10 18 25 16;11 19 26 17;12 20 27 18;13 21 28 19;20 29 29 20;

14 23 31 23;15 24 32 22;16 25 33 23;17 26 34 24;18 27 35 25;19 28 36 26;20 29 37 27;21 30 38 28;29 39 39 29]

%Initialize solution
u=60*ones(1,44);
u(40)=100; u(41)=100; u(42)=100; u(43)=100; u(44)=100;

%Solve the linear system with SOR (Gauss-Seidel=> w=1)

flops(0);

while maxer > .00001
count=count+1;
maxer=0;

for i=1:30
%r=(u(M(i,1))+u(M(i,2))+u(M(i,3))+u(M(i,4)))/4;
r=(1-w)*u(i)+w*(u(M(i,1))+u(M(i,2))+u(M(i,3))+u(M(i,4)))/4;
maxer=max(maxer,abs(r-u(i)));
u(i)=r;

end
for i=1:8
u(30+i)=2*h*(0.2)*(20-u(i+21))+u(i+13);
end
u(39)=2*h*(0.2)*(20-u(30))+u(29);
end

count,flops
u;

%plot the result
x=0:.25:2;
y=0:.25:1;
z=[100,100,100,100,100,u(5),u(11),u(18),u(26);u(1:6),u(12),u(19),u(27);u(7:13),u(20),u(28);
u(14:21),u(29);u(22:30)]
hold off
surf(x,y,z)
view(150,20)
colorbar
rotate3d
view(160,20)

Using SOR on a finer grid

clear
%Irregular Geometry Finite Diff Elliptic PDE Example
h=0.125; n=108; maxer=1; count=0;w=1.6

%Molecule vertex list

M=[126 2 11 2;127 3 12 1;128 4 13 2;129 5 14 3;130 6 15 4;131 7 16 5;132 8 17 6;133 9 18 7;134 10 19 8;9 20 20 9;

1 12 22 12;2 13 23 11;3 14 24 12;4 15 25 13;5 16 26 14;6 17 27 15;7 18 28 16;8 19 29 17;9 20 30 18;10 21 31 19;20 32 32 20;

11 23 34 23;12 24 35 22;13 25 36 23;14 26 37 24;15 27 38 25;16 28 39 26;17 29 40 27;18 30 41 28;19 31 42 29;20 32 43 30;21 33 44 31;32 45 45 32;

22 35 47 35;23 36 48 34;24 37 49 35;25 38 50 36;26 39 51 37;27 40 52 38;28 41 53 39;29 42 54 40;30 43 55 41;31 44 56 42;32 45 57 43;33 46 58 44;45 59 59 45;

34 48 61 48;35 49 62 47;36 50 63 48;37 51 64 49;38 52 65 50;39 53 66 51;40 54 67 52;41 55 68 53;42 56 69 54;43 57 70 55;44 58 71 56;45 59 72 57;46 60 73 58;59 74 74 59;

47 62 76 62;48 63 77 61;49 64 78 62;50 65 79 63;51 66 80 64;52 67 81 65;53 68 82 66;54 69 83 67;55 70 84 68;56 71 85 69;57 72 86 70;58 73 87 71;59 74 88 72;60 75 89 73;74 90 90 74;

61 77 92 77;62 78 93 76;63 79 94 77;64 80 95 78;65 81 96 79;66 82 97 80;67 83 98 81;68 84 99 82;69 85 100 83;70 86 101 84;71 87 102 85;72 88 103 86;73 89 104 87;74 90 105 88;75 91 106 89;90 107 107 90;

76 93 109 93;77 94 110 92;78 95 111 93;79 96 112 94;80 97 113 95;81 98 114 96;82 99 115 97;83 100 116 98;84 101 117 99;85 102 118 100;86 103 119 101;87 104 120 102;88 105 121 103;89 106 122 104;90 107 123 105;91 108 124 106;107 125 125 107];

%Initialize solution

u=80*ones(1,134);
u(126:134)=100;

%Solve the linear system with SOR (Gauss-Seidel=> w=1)

flops(0);
while maxer > .00001
count=count+1;
maxer=0;

for i=1:108
r=(1-w)*u(i)+w*(u(M(i,1))+u(M(i,2))+u(M(i,3))+u(M(i,4)))/4;
maxer=max(maxer,abs(r-u(i)));
u(i)=r;

end

%[count, maxer]

for i=92:107
u(i+17)=2*h*(0.2)*(20-u(i))+u(i-16);
end
u(125)=2*h*(0.2)*(20-u(108))+u(107);
end

count,flops
u;

%plot the result

x=0:.125:2;
y=0:.125:1;

z=[100,100,100,100,100,100,100,100,100,u(9),u(19),u(30),u(42),u(55),u(69),u(84),u(100);
u(1:10),u(20),u(31),u(43),u(56),u(70),u(85),u(101);
u(11:21),u(32),u(44),u(57),u(71),u(86),u(102);
u(22:33),u(45),u(58),u(72),u(87),u(103);
u(34:46),u(59),u(73),u(88),u(104);
u(47:60),u(74),u(89),u(105);
u(61:75),u(90),u(106);
u(76:91),u(107);
u(92:108)]
hold off
surf(x,y,z)
view(150,20)
colorbar
rotate3d

 TOP BACK

HEAT EQUATION

General 1-dimensional Heat Equation Solver using three different methods for comparison

```%Heat Equation Solvers
clear, clf reset
%input parameters
a=0, b=1, k=1, Tfinal=.1, n=20
Const=5 %C=k*DeltaT/DeltaX^2
DeltaX=(b-a)/n
DeltaT=(Const*DeltaX^2)/k
% Initial condition
Uinitial=inline('exp(-((x-(b-a)/2)*100).^2)','x','a','b')
%Uinitial=inline('0*x+1','x','a','b')
%Boundary Conditions
Uleft=inline('0*t','t')
Uright=inline('0*t','t')
%Uleft=inline('0*t+1','t')
%Uright=inline('cos(8*t)','t')
%begin algorithm
t=0; count=1; X=linspace(a,b,n+1);T=[t];
Uout(:,1)=Uinitial(X,a,b)';
while 0
%Crank-Nicolson Matrices
A=diag((-Const/2)*ones(n-2,1),-1)+diag((1+Const)*ones(n-1,1))+diag((-Const/2)*ones(n-2,1),1);
B=diag((Const/2)*ones(n-2,1),-1)+diag((1-Const)*ones(n-1,1))+diag((Const/2)*ones(n-2,1),1);
C=inv(A); D=C*B; MaxEigenvalue=norm(eig(D),'inf')
while T(count)<Tfinal & count<1000
%Crank-Nicolson
U=C*(B*Uout(2:n,count) + (Const/2)*[Uleft(T(count))+Uleft(T(count)+DeltaT);zeros(n-3,1);...
Uright(T(count))+Uright(T(count)+DeltaT)]);
%update solution
T=[T,T(count)+DeltaT]; count=count+1;
Uout(:,count)=[Uleft(T(count));U;Uright(T(count))];
end
end
while 0
%Leap Frog Matrices
A=diag((-Const)*ones(n-2,1),-1)+diag((1+2*Const)*ones(n-1,1))+diag((-Const)*ones(n-2,1),1);
B=2*Const*(diag(ones(n-2,1),-1)+diag((-2)*ones(n-1,1))+diag(ones(n-2,1),1));
%Start-up
U=A\(Uout(2:n,1)+Const*[Uleft(DeltaT);zeros(n-3,1);Uright(DeltaT)]);
T=[T,T(count)+DeltaT]; count=count+1;
Uout(:,count)=[Uleft(T(count));U;Uright(T(count))];
while T(count)<Tfinal & count<1000
%Leap Frog
U=B*Uout(2:n,count)+Uout(2:n,count-1)+(2*Const)*[Uleft(T(count));zeros(n-3,1);Uright(T(count))];
%update solution
T=[T,T(count)+DeltaT]; count=count+1;
Uout(:,count)=[Uleft(T(count));U;Uright(T(count))];
end
end
%while 0
%BDF Matrices
A=diag((-Const)*ones(n-2,1),-1)+diag((1+2*Const)*ones(n-1,1))+diag((-Const)*ones(n-2,1),1);
B=(-2/3*Const)*diag(ones(n-2,1),-1)+diag((1+4*Const/3)*ones(n-1,1))+(-2/3*Const)*diag(ones(n-2,1),1);
C=inv(B);
%Start-up
U=A\(Uout(2:n,1)+Const*[Uleft(DeltaT);zeros(n-3,1);Uright(DeltaT)]);
T=[T,T(count)+DeltaT]; count=count+1;
Uout(:,count)=[Uleft(T(count));U;Uright(T(count))];
while T(count)<Tfinal & count<1000
%BDF Method
U=C*(4/3*Uout(2:n,count)-1/3*Uout(2:n,count-1)+(2*Const/3)*[Uleft(T(count)+DeltaT);zeros(n-3,1);Uright(T(count)+DeltaT)]);
%update solution
T=[T,T(count)+DeltaT]; count=count+1;
Uout(:,count)=[Uleft(T(count));U;Uright(T(count))];
end
%end
count
surf(X,T,Uout'),xlabel('x-axis'),ylabel('t-axis')
colormap([.9,.9,.9])
rotate3d on

 TOP BACK

WAVE EQUATION

```%Wave Equation Solvers
clear, clf reset
%input parameters
a=0, b=1, c=1, Tfinal=2, n=200
DeltaX=(b-a)/n
DeltaT=.005
Const=(c*DeltaT/DeltaX)^2
% Initial condition
Uinitial=inline('exp(-((x-(b+a)/2)*100).^2)','x','a','b')
dUinitial=inline('0*x','x','a','b')
%Boundary Conditions
Uleft=inline('0*t','t')
Uright=inline('0*t','t')
%begin algorithm
t=0; count=1; X=linspace(a,b,n+1);T=[t];
Uout(:,1)=Uinitial(X,a,b)'; Temp=dUinitial(X,a,b)';dU(:,1)=Temp(2:n);
%Leap Frog Matrices
A=diag((Const)*ones(n-2,1),-1)+diag((2-2*Const)*ones(n-1,1))+diag((Const)*ones(n-2,1),1);
%Start-up
U=Uout(2:n,1)+DeltaT*dU(:,1);
T=[T,T(count)+DeltaT]; count=count+1;
Uout(:,count)=[Uleft(T(count));U;Uright(T(count))];
while T(count)<Tfinal & count<1000
%Leap Frog
U=A*Uout(2:n,count)-Uout(2:n,count-1)+(Const)*[Uleft(T(count));zeros(n-3,1);Uright(T(count))];
%update solution
T=[T,T(count)+DeltaT]; count=count+1;
Uout(:,count)=[Uleft(T(count));U;Uright(T(count))];
end
count
for j=1:length(T)
figure(1)
plot(X,Uout(:,j)','b')
axis([a,b,-1,1])
text(.3,-.6,['t= ',num2str(T(j),3)])
M(j) = getframe;
end
movie(M,1)
figure(2)
surf(X,T,Uout')
%colormap([.9,.9,.9])
rotate3d on

 TOP BACK

LINEAR SUBDOMAIN

Code for Example 3.1 followed by Code for Example 3.2

%Solves Example 3.1 Using Subdomain Method

clear all, clf reset

% n=#of nodes

n=3, A=[];

for i=1:6

deltax=1/(n-1)

% Build the n by n Tridiagonal System

LowerD= (-1/2+deltax/8)*ones(1,n-1);
MainD=[1, (3/4)*deltax*ones(1,n-2), 3/8*deltax + 1/2];
UpperD=[0, (1/2+deltax/8)*ones(1,n-2)];
B=[2, deltax*ones(1,n-2), deltax/2];
u=tridiag(LowerD, MainD, UpperD, B);
x=linspace(0,1,n);
plot(x,u)
xlabel('x-axis'), ylabel('u-axis'), title('Example 3.1')
grid on, hold on
A=[A;n,u(n)];
n=(n-1)*2+1

end

hold off
[nn,mm]=size(A);
diff(1)=0;ratio(1)=0; ratio(2)=0;

for i=1:nn-1,

diff(i+1)=A(i+1,2)-A(i,2);

end

for i=2:nn-1,

ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,2) + diff(i)/(1-ratio(i+1));

end

z=[A(:,1)';A(:,2)'; diff; ratio; xacc];
fprintf(1,'\n\n n y diff ratio y acc\n\n')
fprintf(1,'%3.0f %16.14f %13.10f %13.10f %16.14f\n',z)

Code for Example 3.2

%Solves Example 3.2 Using Subdomain Method

clear all, clf reset

% n=#of nodes

n=3, A=[];

for i=1:6

deltax=1/(n-1)

% Build the n by n Tridiagonal System

LowerD= [(1/deltax+deltax/8)*ones(1,n-2), 0];
MainD=[1, (-2/deltax+(3/4)*deltax)*ones(1,n-2), 1];
UpperD=[0, (1/deltax+deltax/8)*ones(1,n-2)];
B=[1, zeros(1,n-1)];
u=tridiag(LowerD, MainD, UpperD, B);
x=linspace(0,1,n);
plot(x,u)
xlabel('x-axis'), ylabel('u-axis'), title('Example 3.1')
grid on, hold on
A=[A;n,u((n+1)/2)];
n=(n-1)*2+1

end

hold off
[nn,mm]=size(A);
diff(1)=0;ratio(1)=0; ratio(2)=0;

for i=1:nn-1,

diff(i+1)=A(i+1,2)-A(i,2);

end

for i=2:nn-1,

ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=A(i-1,2) + diff(i)/(1-ratio(i+1));

end

z=[A(:,1)';A(:,2)'; diff; ratio; xacc];
fprintf(1,'\n\n n y diff ratio y acc\n\n')
fprintf(1,'%3.0f %16.14f %13.10f %13.10f %16.14f\n',z)

 TOP BACK

COLLOCATION
Code for the equation du/dx + u =1, u(0) = 2 using Piecewise Linear Basis Functions

clear
hold off
a=0, b=1, n=6, u0=2
Dx=(b-a)/(n-1)
d=[1,(1/2+1/Dx)*ones(1,n-1)]
l=(1/2-1/Dx)*ones(1,n-1)

%build the matrix A
A=diag(d)+diag(l,-1)
B=[u0, ones(1,n-1)]
y=A\B'
x=linspace(a,b,n)
plot(x,y)
grid on

Code for Piecewise Cubic Hermite Basis Functions
Note that there are two function subprograms required to run the Main Program

MAIN PROGRAM
Parameters set for the equation d2u/dx2 + u = 0, u(0) = 1, u(1) = 0

%collocation

clear, clf reset
a=0, b=1, nn=3, AA=[];

for kk=1:6

x=linspace(a,b,nn);
n=length(x);
A=zeros(2*n);
B=zeros(2*n,1);

%boundary conditions

A(1,1)=1; B(1)=1;
A(2*n,2*n-1)=1; B(2*n)=0;

%for each element

for i=1:n-1,

M=hcoll(x(i),x(i+1));
A(2*i:2*i+1,2*i-1:2*i+2)=M;

end

u = A\B;

for i=1:n,

v(i)=u(2*i-1); dv(i)=u(2*i);

end

[x',v'];
hold on
pchpltv2(x,v,dv)
grid on, xlabel('r axis'), ylabel('temperature axis')
AA=[AA;n,v((n+1)/2)];
nn=(nn-1)*2+1

end

hold off
[nn,mm]=size(AA);
diff(1)=0;ratio(1)=0; ratio(2)=0;

for i=1:nn-1,

diff(i+1)=AA(i+1,2)-AA(i,2);

end

for i=2:nn-1,

ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=AA(i-1,2) + diff(i)/(1-ratio(i+1));

end

z=[AA(:,1)';AA(:,2)'; diff; ratio; xacc];
fprintf(1,'\n\n n y diff ratio y acc\n\n')
fprintf(1,'%3.0f %16.14f %13.10f %13.10f %16.14f\n',z)

ELEMENT MATRIX BUILDER CODE

function M = hcoll(p,q)

%computes the matrix entries for the (p,q) element

a=(p+q)/2; r=(q-p)/2;
v=[1, 1/r, 1/(r)^2]';
R=[v, r.*v, v, r.*v];

H1=[.88490017946 .262891711532 .11509982054 -.0704416218017
-.5 .288675134595 .5 -.288675134595
-.86602540378444 -1.36602540378444 .86602540378444 -.36602540378444];

H2=[.11509982054 .0704416218017 .88490017946 -.262891711532
-.5 -.288675134595 .5 .288675134595
.86602540378444 .36602540378444 -.86602540378444 1.36602540378444];

%p1 and p2 are the points in (p,q) that correspond to +-1/sqrt(3)

p1=a-r/sqrt(3); p2=a+r/sqrt(3);

%vector [A,B,C] elements are the coef of the linear ode

%A function, B first deriv, C second deriv

M(1,:)=[1, 0, 1]*(H1.*R);
M(2,:)=[1, 0, 1]*(H2.*R);

PIECEWISE CUBIC HERMITE GRAPHING SUBPROGRAM

function pchpltv2(x,u,du)

n=length(x);
y=[]; t=[];

for i=1:n-1,

Dx=x(i+1)-x(i); Du=u(i+1)-u(i);FD=Du/Dx;
A=(FD-du(i))/Dx; SD=(du(i+1)-FD)/Dx; B=(SD-A)/Dx;
p=linspace(x(i),x(i+1),16);
s=u(i)+(p-x(i)).*(du(i)+(p-x(i)).*(A+(p-x(i+1)).*B));
y=[y,s]; t=[t,p];
[u(i),du(i), A, B];

end

plot(t,y,'b')

 TOP BACK

GALERKIN

%PLGalerkin for the Linear Second Order Equation
%cC d2u/dx2 + cB du/dx +cA u = f on the interval [a,b]

clear, clf reset
a=0, b=1, nn=3, AA=[]; cC=1,cB=0,cA=1

%Master element constants

I=[2/3 1/3;1/3 2/3]; DI = [-1/2 1/2; -1/2 1/2]; DD = [1/2 -1/2; -1/2 1/2];

for kk=1:6, %for 6 grid refinements

x=linspace(a,b,nn);
n=length(x);

%Step #1 set A matrix and B vector to all zeros

A = zeros(n); B = diag(A); Z=B;

%Step #2 build A matrix from the 2x2 blocks

for i=1:n-1, %for each element

%compute radius and midpoint of ith element

r = (x(i+1)-x(i))/2; ave=x(i)+r;

%build and load the 2x2 block

iI=[i,i+1];
A(iI,iI)=A(iI,iI)+cA*I.*r+cB*DI-cC*DD./r;

end

%Step #3 build the B vector

%B=2*r*ones(n,1); B(1)=r;B(n)=r;%Example 1

%B=r*[(x(1)+x(2))/2, (x(1:n-2)+2*x(2:n-1)+x(3:n))/2, (x(n-1)+x(n))/2]';%Ex 4

%Step #4 correct for the boundary conditions

A(1,:)=Z'; %puts zeros in the first row

A(1,1) = 1; B(1) = 1;
A(n,:)=Z'; A(n,n)=1; B(n)=0; %Ex 2

%A(n,n)=A(n,n)-1/2;B(n)=B(n)-3;%Ex 3

%Step #5 solve the system

u=A\B;
hold on
plot(x,u,'b')
grid on, xlabel('x axis'),ylabel('u axis')
AA=[AA;n,u((n+1)/2)];
nn=(nn-1)*2+1;

end

hold off
[nn,mm]=size(AA);
diff(1)=0;ratio(1)=0; ratio(2)=0;xacc(1)=0;xacc(2)=0;

for i=1:nn-1,

diff(i+1)=AA(i+1,2)-AA(i,2);

end

for i=2:nn-1,

ratio(i+1)=diff(i+1)/diff(i);
xacc(i+1)=AA(i-1,2) + diff(i)/(1-ratio(i+1));

end

z=[AA(:,1)';AA(:,2)'; diff; ratio; xacc];
fprintf(1,'\n\n n y diff ratio y acc\n\n')
fprintf(1,'%3.0f %16.14f %13.10f %13.10f %16.14f\n',z)

 TOP BACK

TRIANGULAR GALERKIN

%TRIANGULARGalerkin

% coordinates of the nodes

clear, clf reset
p=pi/12;
x=[1 2 cos(2*p) 3/2 2*cos(p) 2*cos(2*p) 1 2*cos(3*p) cos(4*p) 1/2 2*cos(4*p) 2*cos(5*p) 0 0];
y=[0 0 sin(2*p) 1/2 2*sin(p) 2*sin(2*p) 1 2*sin(3*p) sin(4*p) 3/2 2*sin(4*p) 2*sin(5*p) 1 2];
nn=length(x);

%incidence matrix

M=[1 2 4;1 4 3;2 5 4;3 4 7;4 6 7;4 5 6;3 7 9;7 6 8;9 7 10;7 11 10;7 8 11;9 10 13;10 11 12;13 10 14;10 12 14];
[n,r]=size(M);

%Step #1 set A matrix and B vector to all zeros

A = zeros(nn); B = diag(A);

%Step #2 build A matrix for each element

for i=1:n,

xi=[x(M(i,:))]; yi=[y(M(i,:))];
a=det([1 1 1;xi;yi])/2;
Mi=[M(i,:) M(i,:)]; xi=[xi xi xi]; yi=[yi yi yi];

for j=1:3,

for k=0:2,

l=j+k;
A(Mi(j),Mi(l))=A(Mi(j),Mi(l))-...
((yi(j+1)-yi(j+2))*(yi(l+1)-yi(l+2))+(xi(j+2)-xi(j+1))*(xi(l+2)-xi(l+1)))/(4*a);

end

end

end

%Step #3 build the B vector

%Step #4 correct for the boundary conditions

%type 1 boundary conditions (replace that row with the condition equation)

Z=zeros(1,nn); B=Z';
A(1,:)=Z; A(1,1)=1; B(1)=2;A(2,:)=Z; A(2,2)=1; B(2)=3;
A(3,:)=Z; A(3,3)=1; B(3)=2;A(9,:)=Z; A(9,9)=1; B(9)=2;
A(13,:)=Z; A(13,13)=1; B(13)=2;A(14,:)=Z; A(14,14)=1; B(14)=3;

%type 2 boundary conditions (adjust the B vector)

L=sqrt((x(2)-x(5))^2+(y(2)-y(5))^2);
B(5)=B(5)-L; B(6)=B(6)-L;B(8)=B(8)-L; B(11)=B(11)-L;B(12)=B(12)-L;

%type 3 boundary conditions (adjust both the B vector and A matrix)

%A(1,1)=A(1,1)-.125; A(1,2)=A(1,2)-.125; B(1)=B(1)-.25;%Example #4
%A(2,1)=A(2,1)-.125; A(2,2)=A(2,2)-.125; B(2)=B(2)-.25;

%Step #5 solve the system

u=A\B;
[x;y;u'], spy(A), pause
trisurf(M,x,y,zeros(1,nn),'EdgeColor','g','LineWidth',2)
hold on
trisurf(M,x,y,u','FaceColor','interp','EdgeColor','b','LineWidth',1)
xlabel('x axis'),ylabel('y axis'),zlabel('u axis'), colorbar
grid on, rotate3d

 TOP BACK

Triangular Grid Refinement

```%Grid Refinement Program
% coordinates of the nodes
clear, clf reset
p=pi/12;
x=[1 2 cos(2*p) 3/2 2*cos(p) 2*cos(2*p) 1 2*cos(3*p) cos(4*p) 1/2 2*cos(4*p) 2*cos(5*p) 0 0];
y=[0 0 sin(2*p) 1/2 2*sin(p) 2*sin(2*p) 1 2*sin(3*p) sin(4*p) 3/2 2*sin(4*p) 2*sin(5*p) 1 2];
nn=length(x);
%incidence matrix
M=[1 2 4;1 4 3;2 5 4;3 4 7;4 6 7;4 5 6;3 7 9;7 6 8;9 7 10;7 11 10;7 8 11;9 10 13;10 11 12;13 10 14;10 12 14];
%boundary point vectors
B1=[1 3 9 13];B2=[2 5 6 8 11 12 14];B3=[1 2];B4=[13 14];
[n,r]=size(M)
figure(1)
trisurf(M,x,y,zeros(1,nn),'FaceColor','none','EdgeColor','g','LineWidth',2)
for i=1:length(x)
hold on
text(x(i),y(i),num2str(i),'fontweight','bold')
end
%Begin grid refinement
Mnew=[];
for i=1:n
MM=[M(i,:),M(i,:)];
for j=1:3
midj=[(x(MM(j))+x(MM(j+1)))/2,(y(MM(j))+y(MM(j+1)))/2];
if ~(isempty(find(MM(j)==B1)))&~(isempty(find(MM(j+1)==B1))),%Boundary #1
rcorr=1/sqrt(midj*midj'); midj=rcorr*midj;
x=[x,midj(1)]; y=[y,midj(2)];
ind(j)=length(x); B1=[B1,ind(j)];
elseif ~(isempty(find(MM(j)==B2)))&~(isempty(find(MM(j+1)==B2))),%Boundary #2
rcorr=1/sqrt(midj*midj'); midj=2*rcorr*midj;
x=[x,midj(1)]; y=[y,midj(2)];
ind(j)=length(x); B2=[B2,ind(j)];
elseif ~(isempty(find(MM(j)==B3)))&~(isempty(find(MM(j+1)==B3))),%Boundary #3
x=[x,midj(1)]; y=[y,midj(2)];
ind(j)=length(x); B3=[B3,ind(j)];
elseif ~(isempty(find(MM(j)==B4)))&~(isempty(find(MM(j+1)==B4))),%Boundary #4
x=[x,midj(1)]; y=[y,midj(2)];
ind(j)=length(x); B4=[B4,ind(j)];
else
ix=find(x==midj(1)*ones(1,length(x))&y==midj(2)*ones(1,length(y))),%repeated points?
if isempty(ix)
x=[x,midj(1)];
y=[y,midj(2)];
ind(j)=length(x);
else
ind(j)=ix;
end
end
end
Mnew=[Mnew;MM(1),ind(1),ind(3);MM(2),ind(2),ind(1);MM(3),ind(3),ind(2);ind(1),ind(2),ind(3)];
end
%put boundary points in order
[So,II]=sort(x(B1));B1=B1(II);
[So,II]=sort(x(B2));B2=B2(II);
[So,II]=sort(x(B3));B3=B3(II);
[So,II]=sort(y(B4));B4=B4(II);
M=Mnew;
%End of grid refinement
nn=length(x);
figure(2)
trisurf(M,x,y,zeros(1,nn),'FaceColor','none','EdgeColor','g','LineWidth',2)
[x;y],B1,B2,B3,B4
for i=1:length(x)
hold on
text(x(i),y(i),num2str(i),'fontweight','bold')
end```

 TOP BACK

Revised 12/06/04
By Don Short, dshort@sciences.sdsu.edu