SDSU_signature-300.gif (1279 bytes)

MATH 542 NUMERICAL SOLUTIONS OF DIFFERENTIAL EQUATIONS


INTRODUCTION
INITIAL VALUE PROBLEMS(ODE)
FINITE DIFFERENCE METHODS
FINITE ELEMENT METHOD
HOME

GRADIENT FIELD WITH SOLUTIONS
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(Heading)
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(Heading)
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(Heading)

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)
shading interp, view(150,30),colorbar
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)
shading interp
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)
shading interp
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

ADAMS METHODS

2nd Order Adams Method

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


4th Order Adams Method

function [u,v] = adam4(a,b,h,ya)
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)
shading interp
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)
shading interp
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)
shading interp
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)
shading interp
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
%shading interp

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
shading interp

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