Matlab implementation of the examples presented in the lecture notes of Numerieke Wiskunde voor technici

• Chapter 1 Beginwaardenproblemen
• Chapter 2 Het oplossen van stelsels lineaire vergelijkingen

1. Beginwaardenproblemen

1. Example 1.4 page 17
2. Example 1.7 page 32
3. Example 1.8 page 37
4. Example 1.10 page 46

 
%
% This matlab program illustrates example 1.4 from Numerieke Wiskunde
% voor Technici. The model describes the flow of a fluid from a
% reservoir.
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 18-11-1997

clf
hold off
clear

n = 10;
h = 5/n;

for i = 1 :n+1
   x(i) = (i-1)*h;
end

y(1) = 0;
a = 1;
f = 1;

for i = 1 :n+1
   y(i) = sqrt(f/a)*tanh(x(i)*sqrt(a*f));
end

yex = y;
hp = plot(x,y,'k');
set(hp, 'LineWidth',2);
hold on

% Euler

for i = 1 :n
   y(i+1) = y(i) + h*(f-a*(y(i)*y(i)));
end
hp = plot(x,y,'*r');
set(hp, 'LineWidth',2);
max(abs(y-yex))

% Heun

for i = 1 :n
   y(i+1) = y(i) + h*(f-a*(y(i)*y(i)));
   y(i+1) = y(i) + h*(f-a*(y(i)*y(i))+f-a*(y(i+1)*y(i+1)))/2;
end
max(abs(y-yex))
hp = plot(x,y,'ob');
set(hp, 'LineWidth',2);

% Runge Kutta

for i = 1 :n
   k1 = h*(f-a*(y(i)*y(i)));
   k2 = h*(f-a*((y(i)+k1/2)^2));
   k3 = h*(f-a*((y(i)+k2/2)^2));
   k4 = h*(f-a*((y(i)+k3)^2));
   y(i+1) = y(i) + (k1+2*k2+2*k3+k4)/6;
end
max(abs(y-yex))
hp = plot(x,y,'+g');
set(hp, 'LineWidth',2);

xlabel('t','FontSize', 16);
ylabel('Debiet Q','FontSize', 16)
title(['Oplossing voor ' num2str(n) ' punten'],'FontSize', 16)
x=[3.5 4];
y=[0.1 0.1];
hp = plot(x,y,'k');
set(hp, 'LineWidth',2);
text(4.1,0.1,'Exact','FontSize', 16);
y=[0.2 0.2];
hp = plot(x,y,'*r');
set(hp, 'LineWidth',2);
text(4.1,0.2,'Euler','FontSize', 16);
y=[0.3 0.3];
hp = plot(x,y,'ob');
set(hp, 'LineWidth',2);
text(4.1,0.3,'Heun','FontSize', 16);
y=[0.4 0.4];
hp = plot(x,y,'+g');
set(hp, 'LineWidth',2);
text(4.1,0.4,'Runge Kutta','FontSize', 16);

 

 
%
% This matlab program illustrates example 1.7 from Numerieke Wiskunde
% voor Technici. The model describes the flow of a fluid from a
% reservoir. The (in)stability of the Euler method is investigated.
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 21-11-1997

clf
hold off
clear

n = 12;
h = 6/n;

for i = 1 :n+1
   x(i) = (i-1)*h;
end

y(1) = 0;
x(1) = 0;
a = 1;
f = 1;

for i = 1 :n+1
   y(i) = sqrt(f/a)*tanh(x(i)*sqrt(a*f));
end

yex = y;
hp = plot(x,y,'k');
set(hp, 'LineWidth',2);
axis([0 6 0 1.5]);
hold on

clear x y
y(1) = 0;
x(1) = 0;
for i = 1 :n
   y(i+1) = y(i) + h*(f-a*(y(i)*y(i)));
   x(i+1) = i*h;
end

hp = plot(x,y,'g');
set(hp, 'LineWidth',2);

clear x y
y(1) = 0;
x(1) = 0;
n = 6;
h = 6/n;
for i = 1 :n
   y(i+1) = y(i) + h*(f-a*(y(i)*y(i)));
   x(i+1) = i*h;
end

hp = plot(x,y,'b');
set(hp, 'LineWidth',2);

clear x y
y(1) = 0;
x(1) = 0;
n = 3;
h = 6/n;
for i = 1 :n
   y(i+1) = y(i) + h*(f-a*(y(i)*y(i)));
   x(i+1) = i*h;
end

hp = plot(x,y,'r');
set(hp, 'LineWidth',2);

clear x y
xlabel('t','FontSize', 16);
ylabel('Debiet Q','FontSize', 16)
title(['Oplossing voor Debiet probleem'],'FontSize', 16)
x=[3.5 4.2];
y=[0.1 0.1];
hp = plot(x,y,'k');
set(hp, 'LineWidth',2);
text(5.0,0.1,'Exact','FontSize', 16);
y=[0.2 0.2];
hp = plot(x,y,'g');
set(hp, 'LineWidth',2);
text(5.0,0.2,'Euler h = 0.5','FontSize', 16);
y=[0.3 0.3];
hp = plot(x,y,'b');
set(hp, 'LineWidth',2);
text(5.0,0.3,'Euler  h = 1','FontSize', 16);
y=[0.4 0.4];
hp = plot(x,y,'r');
set(hp, 'LineWidth',2);
text(5.0,0.4,'Euler  h = 2','FontSize', 16);

 

 

%
% This matlab program illustrates example 1.8 from Numerieke Wiskunde
% voor Technici. The model describes a predator-prey problem. The
% resulting system of differential equations is non-linear.
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 25-11-1997

clear u v t;
hold off
clf
n = 250;
us = 100;
vs = 2;
b=1000;
a=1/40000;
k = 6000/n;
u(1) = us;
v(1) = vs;
t(1) = 0;
for j = 1 : n;
   u(j+1) = u(j)+k*(a*(b-u(j))-(v(j)/200))*u(j);
   v(j+1) = v(j)+k*(0.000023*u(j)-0.0006)*v(j);
   t(j+1) = j*k;
end;
hp = plot(t,u,'xr');
set(hp, 'LineWidth',2);
hold on
x = [4000 4500];
y = [300 300];
hp = plot(x,y,'xr');
set(hp, 'LineWidth',2);
text(5000,300,'h = 24','FontSize', 16);

title('Roofdier-prooi model opgelost met Euler','FontSize', 16)
xlabel('tijd','FontSize', 16)
ylabel('Prooidieren','FontSize', 16)

clear u v t;
n = 500;
k = 6000/n;
u(1) = us;
v(1) = vs;
t(1) = 0;
for j = 1 : n;
   u(j+1) = u(j)+k*(a*(b-u(j))-(v(j)/200))*u(j);
   v(j+1) = v(j)+k*(0.000023*u(j)-0.0006)*v(j);
   t(j+1) = j*k;
end;
hp = plot(t,u,'g');
set(hp, 'LineWidth',2);
x = [4000 4500];
y = [250 250];
hp = plot(x,y,'g');
set(hp, 'LineWidth',2);
text(5000,250,'h = 12','FontSize', 16);

clear u v t;
n = 1000;
k = 6000/n;
u(1) = us;
v(1) = vs;
t(1) = 0;
for j = 1 : n;
   u(j+1) = u(j)+k*(a*(b-u(j))-(v(j)/200))*u(j);
   v(j+1) = v(j)+k*(0.000023*u(j)-0.0006)*v(j);
   t(j+1) = j*k;
end;
hp = plot(t,u,'b');
set(hp, 'LineWidth',2);
x = [4000 4500];
y = [200 200];
hp = plot(x,y,'b');
set(hp, 'LineWidth',2);
text(5000,200,'h = 12','FontSize', 16);

 

 

%
% This matlab program illustrates example 1.10 from Numerieke Wiskunde
% voor Technici. We solve the mathematical pendulum with the Euler
% forward and the Runge Kutta method. Since the system of linearized
% equations has eigenvalues on the imaginary axis, only the solution
% obtained with Runge Kutta is stable.
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 25-11-1997

clear x y xhulp yhulp t
hold off
clf

% Euler

n = 1000;
dt = 1/10;
x(1) = 1;
y(1) = 0;
t(1) = 0;
for j = 1:n
   t(j+1) = j*dt;
   x(j+1) = x(j)+dt*y(j);
   y(j+1) = y(j)-dt*sin(x(j));
end

hp = plot(t,x,'r');
set(hp, 'LineWidth',2);
hold on

title('De mathematische slinger','FontSize', 16)
xlabel('tijd','FontSize', 16)
ylabel('Hoekuitwijking','FontSize', 16)

% Runge Kutta

x(1) = 1;
y(1) = 0;
for j = 1:n
   k1(1) = dt*y(j);
   k1(2) = -dt*sin(x(j));
   k2(1) = dt*(y(j)+k1(2)/2);
   k2(2) = -dt*sin(x(j)+k1(1)/2);
   k3(1) = dt*(y(j)+k2(2)/2);
   k3(2) = -dt*sin(x(j)+k2(1)/2);
   k4(1) = dt*(y(j)+k3(2));
   k4(2) = -dt*sin(x(j)+k3(1));
   x(j+1) = x(j)+(k1(1)+2*k2(1)+2*k3(1)+k4(1))/6;
   y(j+1) = y(j)+(k1(2)+2*k2(2)+2*k3(2)+k4(2))/6;
end
hp = plot(t,x,'g');
set(hp, 'LineWidth',2);

 

2. Het oplossen van stelsels lineaire vergelijkingen

1. Example 2.6_1 Section 2.6

 
%
% This matlab program illustrates the application given in Section 2.6
% from Numerieke Wiskunde voor Technici. We solve the convection-difussion
% equation, for various values of the velocity v. It appears that for
% small values of v, a central discretization is the best choice, whereas
% for large values of v, upwind discretization leads to the best results.
%
% programmer: C. Vuik, TU Delft, The Netherlands
% e-mail    : c.vuik@math.tudelft.nl
% date      : 16-12-1997


% This is an example of a one-dimensional convection diffusion equation
% 
%        -u'' + v u' = 1  u(0) = 0, u(1) = 0

clear x u een twee a y b
hold off
clf

% change the value of v to see the difference between the central and
% upwind discretization of the problem.

v = 100;
n = 10;
h = 1/n;

% the analytical solution

c2 = -1/(v*(1-exp(v)));
c1 = -c2;

for i = 1: n+1
   x(i) = (i-1)*h;
   u(i) = c1*exp(v*x(i))+c2+x(i)/v;
end
hp = plot(x,u,'g');
set(hp, 'LineWidth',2);
hold on

% central discretization

for i = 1:n-1
   twee(i,i) = 2;
   b(i,1)      = 1;
end
for i = 1:n-2
   twee(i+1,i) = -1;
   twee(i,i+1) = -1;
   een (i+1,i) = -1;
   een (i,i+1) =  1;
end

een = sparse(een);
twee = sparse(twee);

a = twee/(h*h) + v*een/(2*h);

y = a\b;

for i = 1:n-1
   u(i+1)      = y(i);
end

hp = plot(x,u,'*');
set(hp, 'LineWidth',2);

% upwind discretization

clear een

for i = 1:n-1
   een(i,i) = 1;
end
for i = 1:n-2
   een (i+1,i) = -1;
end
een = sparse(een);

a = twee/(h*h) + v*een/(h);

y = a\b;

for i = 1:n-1
   u(i+1)      = y(i);
end

hp = plot(x,u,'ro');
set(hp, 'LineWidth',2);

title ('Convection-diffusion equation','FontSize', 16);
xlabel('x-coordinate','FontSize', 16);
ylabel('Temperature','FontSize', 16);

 

Contact information: