MATH 542 NUMERICAL SOLUTIONS OF DIFFERENTIAL EQUATIONS

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

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

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

clear
global func_f func_g
clf reset

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

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

%enter f and g as a character string

clear

global func_f func_g

clf reset

%set the figure window size and grid size

[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

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)

grid on

%axis labels

xlabel('time axis')

ylabel('angle axis')

zlabel('velocity axis')

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

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)

FIRST ODER PDE EXAMPLE

%solution of a first order PDE

TRAFFIC MODEL W/O SHOCK

%Traffic Model without shock

%calculation of d0 = d(x,0)

for i=1:n,

for j=1:m,

end

end

TRAFFIC MODEL WITH SHOCK

%Traffic Model with shock

%calculation of d0 = d(x,0)

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

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)

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

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)

for i=1:n

end

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);

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

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

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')

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')

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

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)')

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

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

%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

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

%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

%plot the result

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

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

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

% 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

for i=1:nn-1,

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

for i=2:nn-1,

end

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

% 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

for i=1:nn-1,

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

for i=2:nn-1,

end

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

for kk=1:6

%boundary conditions

%for each element

for i=1:n-1,

end

u = A\B;

for i=1:n,

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

end

end

for i=1:nn-1,

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

for i=2:nn-1,

end

ELEMENT MATRIX BUILDER CODE

function M = hcoll(p,q)

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

%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);

%A function, B first deriv, C second deriv

PIECEWISE CUBIC HERMITE GRAPHING SUBPROGRAM

function pchpltv2(x,u,du)

for i=1:n-1,

end

plot(t,y,'b')

GALERKIN

%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

%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

%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(n,n)=A(n,n)-1/2;B(n)=B(n)-3;%Ex 3

%Step #5 solve the system

end

for i=1:nn-1,

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

for i=2:nn-1,

end

TRIANGULAR GALERKIN

% coordinates of the nodes

%incidence matrix

%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,

for j=1:3,

for k=0:2,

end

end

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

%type 2 boundary conditions (adjust the B vector)

%Step #5 solve the system

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

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