Simulation tools

Part: Numerical Methods

Prepared by [Link] Iqteit

2024

1
 SECTION 1: Solving Nonlinear Systems

Newton’s method for solving one equation:

Taylor Series

find x so that f(x) = 0

Then we define a sequence of points {x0, x1, x2, x3, . . .} from the formula:

Newton’s method with number of iteration Newton’s method with certain tolerance for |rn| = |f(xn)|
function x = mynewton (f ,f1 function x = mynewtontol(f,f1,x0,tol)
,x0 ,n) x = x0;
x = x0 ; i =0;
for i = 1: n y = f(x);
x = x - f(x)/f1(x); while abs(y) > tol & i < 1000
fx=f(x); x = x - y/f1(x)
end y = f(x);
xans=x i = i +1;
fxfinal=fx end
end xans=x
fxfinal=y
Example: solve equation end
𝑥𝑒 𝑥 − sin(𝑥) = 10 Example: solve equation
𝑥𝑒 𝑥 − sin(𝑥) = 10
>> f=@(x)x.*exp(x)-sin(x)-10;
>> f1=@(x)exp(x)+x.*exp(x)-cos(x); >> f=@(x)x.*exp(x)-sin(x)-10;
>> x0=2; >> f1=@(x)exp(x)+x.*exp(x)-cos(x);
>> n=5; >> x0=2;
>> x = mynewton (f ,f1 ,x0 ,n) >> tol=0.001;
xans = >> x = mynewtontol (f,f1,x0,tol)
1.8049 xans =
1.8049
fxfinal =
fxfinal =
3.5527e-15
1.8224e-06

2
Secant Methods for solving one equation

The secant method requires two initial approximations x0 and x1.

From x0 and x1 we can determine that the points (x0, y0 = f(x0)) and (x1, y1 = f(x1)). both lie on the graph of
f. Connecting these points gives the (secant) line:

function x = mysecant(f,x0,x1,n)
y0 = f(x0);
y1 = f(x1);
for i = 1:n
x = x1 - (x1-x0)*y1/(y1-y0);
y=f(x);
x0=x1;
y0=y1;
x1=x;
y1=y;
end
xans=x
fxfinal=y
end
Example: solve equation
𝑥𝑒 𝑥 − sin(𝑥) = 10
>> f=@(x)x.*exp(x)-sin(x)-10;
>> x0=2;
>> x1=1;
>> n=5;
>> x = mysecant(f,x0,x1,n)

xans =
1.8049
fxfinal =
5.8895e-05

Exercises
Use Matlab to solve the following equation by Newton's Method and secant method. At least use five
iterations.

𝑥 3 − (𝑥 + 2)−𝑥 − 20 = 0

3
Newton’s method for solving system of equations:
Vector Notation:

Linear Approximation for Vector Functions:

Df(x0) is an n × n matrix whose entries are the various partial derivative of the components of f, evaluated
at x0. Specifically,

Newton’s Method:

For subsequent steps, we have the following process:

Example: use fsolve function to solve this system of equations

𝑥12 + 𝑥22 = 25

𝑥12 𝑥22
+ + 2𝑥1 𝑥2 = 20
16 9
clear all

4
clc
f =@(x)[x(1).^2+x(2).^2-25 ;(x(1).^2/16)+(x(2).^2/9)+2*x(1).*x(2)-20];

x = [-5;-5];
x = fsolve(f,x)
x=
>> syms x1 x2
-4.5875 >> ezplot(x1^2+x2^2-25)
-1.9886 >> hold on
>> ezplot(x1^2/16+x2^2/9+2*x1*x2-20)

Example:
Solve this system of equations by using Newton’s Method:

𝑥12 + 𝑥22 + 𝑥32 = 6

𝑥12 − 𝑥22 + 2𝑥3 = 2

2𝑥12 + 𝑥22 − 𝑥32 = 3

clear all
clc
n =8;
f = @(x)[x(1)^2+x(2)^2+x(3)^2-6 ; x(1)^2-x(2)^2+2*x(3) -2;
2*x(1)^2+x(2)^2-x(3)^2-3];
Df = @(x)[2*x(1) 2*x(2) 2*x(3);2*x(1) -2*x(2) 2; 4*x(1) 2*x(2) -
2*x(3)];
x = [.5;.5;1];
for i = 1: n
Dx = -Df(x)\f(x);
x = x + Dx;
f ( x );
end
xans=x
yfinal=f(x)
xans =
1.0000
1.7321
1.4142
yfinal =
1.0e-15 *
0

5
 0.4441
-0.4441
Example: use Newton’s Method and symbolic variables for solving this system of equations
2𝑥 + 𝑦 + 2𝑧 = 3

𝑥 2 + 𝑦 2 + 2𝑧 2 = 5
𝑒 𝑥𝑦 + 𝑥𝑦 − 𝑥𝑧 = 1
clear all
clc
syms x y z
f1=2*x+y+2*z-3;
f2=x^2 + y^2+ 2*z^2 -5;
f3=exp(x*y)+x*y-x*z-1;
f=[f1;f2;f3];
Df=[diff(f1,x) diff(f1,y) diff(f1,z)
diff(f2,x) diff(f2,y) diff(f2,z)
diff(f3,x) diff(f3,y) diff(f3,z)];
n=5;
xyz = [2.1;-0.27;-0.47];
for i= 1:n
Dxyz=-double(subs(Df,[x,y,z],xyz')\subs(f,[x,y,z],xyz'));
xyz = xyz + Dxyz;
f=double(subs(f,[x;y;z],xyz));
end
xyzfinal=xyz
yfinal=f
xyz =
2.1138
-0.2716
-0.4780

f =
1.0e-03 *
0.0000
0.3354
0.0954

Exercises
Write code of matlab to solve the following nonlinear algebraic system, then find positive solution of x
and y . At least use five iterations.

𝑥 𝑦 − 𝑒 −2𝑥 − 7 = 0
𝑦 𝑥 − 𝑒 −2𝑦+1 − 7 = 0

6
SECTION 2: Functions and Data
Plotting data
In a chemical reaction the concentration level y of the product at time t was measured every half hour.
The following results were found:

>> t1 = 0:.5:2;
>> y1 = [ 0 .19 .26 .29 .31 ];
>> plot ( t1 , y1,'*- ')
>> title('y=f(t)')
>> xlabel('t')
>> ylabel('y')

Polynomial Interpolation

Lagrange interpolation

Example: find the Lagrange polynomial of this data set

>> x=0:3;
>> y=[-5 -6 -1 16];
>> disp([x;y])
0 1 2 3
-5 -6 -1 16

By mathematics:

7
 (𝑥 − 1)(𝑥 − 2)(𝑥 − 3) (𝑥 − 0)(𝑥 − 2)(𝑥 − 3) (𝑥 − 0)(𝑥 − 1)(𝑥 − 3)
𝑃(𝑥) = (−5) + (−6) + (−1)
(0 − 1)(0 − 2)(0 − 3) (1 − 0)(1 − 2)(1 − 3) (2 − 0)(2 − 1)(2 − 3)
(𝑥 − 0)(𝑥 − 1)(𝑥 − 2)
+ (16)
(3 − 0)(3 − 1)(3 − 2)

𝑃(𝑥) = 𝑥 3 − 2𝑥 − 5

By MATLAB:
Lagrange function in matlab Example
function v=polyinterp(x,y,u) >> x=0:3;
n=length(x); >> y=[-5 -6 -1 16];
v=zeros(size(u)); >> u=-0.05:0.01:3.05;
for k=1:n; >> v=polyinterp(x,y,u);
w=ones(size(u)); >> plot(x,y,'o',u,v,'-')
for j=[1:k-1 k+1:n]
>> xlabel('x and u')
w=(u-x(j)./(x(k)-x(j)).*w;
>> ylabel('y and v')
end
v=v+w*y(k);
end >> symx=sym('x');
>>P=polyinterp(x,y,symx)
end P=
(x*(x - 1)*(x - 3))/2 + 5*(x/2 - 1)*(x/3 - 1)*(x - 1) +
(16*x*(x/2 - 1/2)*(x - 2))/3 - 6*x*(x/2 - 3/2)*(x - 2)
%pretty(P)
>> P=simplify(P)
P=
x^3 - 2*x - 5

Power form interpolation

8
Example A Example B
clear all clear all
clc clc
x=0:1:3; t= [0.2 .5 1.0 1.5 6.0];
y=[-5 -6 -1 16] y= [ 0.5 1.19 2.26 .29 6.31];
A=vander(x) yi=[1.19;2.26;0.29;6.31];
c=inv(A)*y' xi=[ 0.5^4 0.5^3 0.5^2 0.5^1;
syms x 1^4 1^3 1^2 1^1;
f= 1.5^4 1.5^3 1.5^2 1.5^1;
c(1)*x.^3+c(2)*x.^2+c(3)*x+c(4) 6^4 6^3 6^2 6^1];
A = a=inv(xi)*yi
syms x
0 0 0 1 f=a(1)*x^4+a(2)*x^3+a(3)*x^2+a(4)*x+0
1 1 1 1 ezplot(f,0,7)
8 4 2 1 hold on
27 9 3 1 plot(t ,y,'*')
c = a =
0.8651
1.0000 -6.4887
0.0000 7.9791
-2.0000 -0.0955
-5.0000 f =
f =
(12847*x^4)/14850 - (32119*x^3)/4950
x^3 +x^2/1125899906842624-2*x - + (43087*x^2)/5400 - (1891*x)/19800
5 >> vpa(f,4)

ans =

0.8651*x^4 - 6.489*x^3 + 7.979*x^2 -

0.09551*x

Cubic spline data interpolation

clear all
close all

9
clc
tt=1
if tt==1
x=0:0.5*pi:2*pi
y = sin(x);
xx=0:0.1:7
yy = spline(x,y,xx);
plot(x,y,'o',xx,yy)
%----------------------------
elseif tt==2
x = 0:.5:6.5;
Y = [sin(x); cos(x)];
xx = 0:.02:6.5;
YY = spline(x,Y,xx);
plot(x,Y(1,:),'o',xx,YY(1,:),'-')
hold on
plot(x,Y(2,:),'o',xx,YY(2,:),':')
%-----------------------------------
else
x = -3:3;
y = [-1 -1 -1 0 1 1 1];
xq1 = -3:.01:3;
p = pchip(x,y,xq1);
s = spline(x,y,xq1);
plot(x,y,'o',xq1,p,'-',xq1,s,'-.')
legend('SamplePoints','pchip','spline','makima','Location','SouthEast'
)
hold off
end
tt = 1 tt = 2 tt =3

Least Squares Fitting

Suppose we have n data points {𝑥𝑖 , 𝑦𝑖 }𝑛𝑖=1 . In the previous sections we learned about interpolation, where
the fitting function f(x) is required to satisfy yi = f(xi) for i = 1, . . . , n. Very often data has a significant
amount of noise, so it makes more sense to only ask for yi ≈ f(xi). The least squares fitting method produces
such an f. The next illustration shows the effects noise can have and how least squares is used.

10
The least squares method requires you to first select what type of f to consider (sometimes called a model)
and what parameters it depends on. For example, you might choose

Given data {𝑥𝑖 , 𝑦𝑖 }𝑛𝑖=1 , we can then determine the error at each point as ei = yi − f(xi). To make f reasonable,
we wish to simultaneously minimize all the errors {e1, e2, . . . , en}.the standard one is to try to minimize
the sum of the squares of the errors. That is, we denote by E the sum of the squares.

Linear least squares

The function f(x) = a2x2 + a1x + a0 is a quadratic in x, but depends linearly on each of a2, a1, and a0.

Matlab’s basic fitting tool uses this process to obtain a n degree polynomial fit whenever the number of
data points is more than n+1
p = polyfit(x,y,n) returns the coefficients for a polynomial p(x) of degree n that is a best fit (in a least-
squares sense) for the data in y. The coefficients in p are in descending powers, and the length of p is n+1
p(x)=p1xn+p2xn−1+...+pnx+pn+1.

11
Example: Drag due to air resistance is proportional to the square of the velocity, i.e. d = kv2.
>> v = 0:1:60;
>> d = .1234* v .^2;
>> dn = d + .4* v .* randn ( size ( v ));
>> plot (v , dn , '* ')

Example : polynomial fitting

clear all
clc
s=1;
if s==1
x = 1:0.1:20;
T = 0.05*x + 1;
Tn = T + .1*randn(size(x));
p= polyfit(x,Tn,4); %p = polyfit(x,y,n),n=4
f = polyval(p,x);
plot(x,Tn ,'.',x,f,'-')

else
v = 0:1:50;
d = .1234* v.^3;
dn = d + .8* v.^2.*randn(size(v));
p= polyfit(v,dn,2); %p = polyfit(x,y,n), n=2
f = polyval(p,v);
plot(v,dn ,'*',v,f,'-')
end
For s=1 and n=4 For s=2 and n=2

p
= p=

0.0000 -0.0003 0.0042 0.0350 1.0004 1.0e+03 *

0.0100 -0.2249 1.0001

f(v) = p(1)*v.^2+p(2)*v+p(3)

12
Example: exponential fitting
x = (0:0.2:5)';
y = 2*exp(-0.2*x) + 0.1*randn(size(x));
f = fit(x,y,'exp1')
plot(f,x,y)
f=

General model Exp1:

f(x) = a*exp(b*x)
Coefficients (with 95% confidence bounds):
a= 2.062 (1.973, 2.152)
b = -0.2011 (-0.2213, -0.1808)

13
SECTION 3: Numerical Integration

Single integral
Method Formulas Program and example
Analytical integral & Exact method by MATLAB
General approximation > syms x
>> f=x*exp(x);
>> int(f)
ans =
exp(x)*(x - 1)
>> int(f,0,5)
>> double(int(f,0,5))
ans =
594.6526
Left-hand point function L = Intg_leftsum(x,y)
n =length(x);
L = 0;
for i = 1:n-1
L = L + y(i)*(x(i+1) - x(i));
end
end
𝟓
Example:Find ∫𝟎 𝒙𝒆𝒙 𝒅𝒙 = 𝒙𝒆𝒙 − 𝒆𝒙 |𝟓𝟎 = 593.6526
>> x=0:0.01:5;
>> y=x.*exp(x);
>> L = Intg_leftsum(x,y)
L=
590.9497
Right-hand point function L = Intg_rightsum(x,y)
n =length(x);
L = 0;
for i = 1:n-1
L = L + y(i+1)*(x(i+1) - x(i));
end
end
𝟓
Example:Find ∫𝟎 𝒙𝒆𝒙 𝒅𝒙 = 𝒙𝒆𝒙 − 𝒆𝒙 |𝟓𝟎 = 593.6526
>> x=0:0.01:5;
>> y=x.*exp(x);
>> L = Intg_rightsum (x,y)
L=
598.3704
Midpoint rule function L = Intg_midpoint(x,f)
n =length(x);
L = 0;
for i = 1:n -1
L = L+f((x(i+1)+x(i))/2)*(x(i+1)-x(i));
end
end

14
 𝟓
Example:Find ∫𝟎 𝒙𝒆𝒙 𝒅𝒙 = 𝒙𝒆𝒙 − 𝒆𝒙 |𝟓𝟎 = 593.6526
>> x=0:0.01:5;
>> f=@(x)x.*exp(x);
>> L = Intg_midpoint(x,f)
L =
594.6489

Trapezoid rule function T = Intg_trap (x,y)

n =length(x);
T = 0;
for i = 1:n -1
T = T+(y(i)+y(i+1))*(x(i+1)-x(i));
end
T = T/2;
End
𝟓
Example:Find ∫𝟎 𝒙𝒆𝒙 𝒅𝒙 = 𝒙𝒆𝒙 − 𝒆𝒙 |𝟓𝟎 = 593.6526
>> x=0:0.01:5;
>> y=x.*exp(x);
>> T=Intg_trap (x,y)
T =
594.6600
Simpson’s Rule function w = mysimpweights(n)
if rem(n ,2) ~= 0
error('n must be even for Simpsons rule ')
end
w = 2* ones (n+1 ,1);
w(1) = 1; w(n+1) = 1;
w (2:2:n)=4;
end
>> w = mysimpweights(6)% weights for
Simpson ’s rule n=6 interval
w =
1
n: number of intervals and even number 4
2
Qusetion: Try to prove Simpson’s Rule 4
2
4
Note: In Matlab there is a built-in command for 1
definite integrals: integral(f, a, b) where the f is a >> w = mysimpweights(5)
function and a and b are the endpoints. The command
uses “adaptive Simpson quadrature”, a form of Error using mysimpweights (line 3)
Simpson’s rule that checks its own accuracy and adjusts n must be even for Simpsons rule.
the grid size where needed. function I=Int_simpsion(f,a,b,n)
w = mysimpweights(n);
>> f = @(x)x.*exp(x); h=(b-a)/n;
>> I = integral (f ,0 ,5) x=[a:h:b];
f(x);
I= I=h*f(x)*w/3
594.6526 end

15
 >> f=@(x)x.*exp(x);
>> a=0; b=5; n=100;
>> I=Int_simpsion(f,a,b,n)

I =
594.6527
line integral clear all
clc
format long
n = 100;
t =linspace(0,2*pi,n+1);
where C is the counter-clockwise curve around the
x =cos(t);
boundary of the region. We can represent such a curve
y =sin(t);
by consecutive points on it, i.e. 𝑥̅ = (x0, x1, x2, . . . , xn−1,
plot(x,y)
xn), and 𝑦̅ = (y0, y1, y2, . . . , yn−1, yn). Since we are
A = 0;
assuming the curve ends where it begins, we require (xn, for i=1:n
yn) = (x0, y0). A = A-(y(i)+y(i+1))*(x(i+1)-x(i));
end
A = A/2
A =
3.139525976465670

Example: use Trapezoid rule

16
clear all
clc
format long
n = 1000;
t =linspace(0 ,1,n +1);
x =1-3*t;
y =3-3*t;
f1=2*x.*y.^2;
f2=cos(pi*y);
plot(x,y)
A1= 0;A2=0;
for i = 1:n
A1 = A1+(f1(i)+f1(i+1))*(x(i+1)-x(i));
A2 = A2+(f2(i)+f2(i+1))*(y(i+1)-y(i));
end
A1=A1/2
A2=A2/2
A = A1+A2
A1 =
-4.500022500000005
A2 =
-4.475586568020162e-16
A=
-4.500022500000005

Error bounds
The trapezoid, midpoint, and Simpson’s rules are all approximations. As with any approximation, before
you can safely use it, you must know how good (or bad) the approximation might be. For these methods
there are formulas that give upper bounds on the error. In other words, the worst case errors for the methods.
These bounds are given by the following statements:

Trapezoid rule

Midpoint rule

Simpson’s rule

In practice K2 and K4 are themselves approximated from the values of f at the evaluation points.
The most important thing in these error bounds is the dependence on h. To emphasize this dependence,
we sometimes use the order notation O( ). The trapezoid and midpoint method errors are O(h2 ), so the
methods have order 2. The Simpson’s rule error is O(h4 ), so it has order 4. If h is just moderately small,
then there is a huge advantage with Simpson’s method.

17
SECTION 4: Double Integrals for Rectangles
Introduction of Plotting Functions of Two Variables

Functions on Rectangular Grids:

Example
clear all
close all
clc
f = @(x,y)(y.^2).*sin(x+y)
[X,Y] =meshgrid(-5:0.1:5 ,-pi:.05*pi:pi);
Z = f(X,Y);
figure (1)
mesh(X,Y,Z)
figure (2)
surf(X,Y,Z)
XX=size(X)
YY=size(Y)
ZZ=size(Z)

f =
function_handle with value:

@(x,y)(y.^2).*sin(x+y)

XX =
41 101
YY =
41 101
ZZ =
41 101

>> X (1:6 ,1:6)

ans =
-5.0000 -4.9000 -4.8000 -4.7000 -4.6000 -4.5000
-5.0000 -4.9000 -4.8000 -4.7000 -4.6000 -4.5000
-5.0000 -4.9000 -4.8000 -4.7000 -4.6000 -4.5000
-5.0000 -4.9000 -4.8000 -4.7000 -4.6000 -4.5000
-5.0000 -4.9000 -4.8000 -4.7000 -4.6000 -4.5000
-5.0000 -4.9000 -4.8000 -4.7000 -4.6000 -4.5000

18
The center point method

If we denote the center points by cij , then the sum becomes

19
 function I=Integralcenterpoint(f,a,b,c,d,n,m)
x=a:(b-a)/n:b;
y=c:(d-c)/n:d;
s=0;
for ii=1:n
for jj=1:m
s=s+f((x(ii)+x(ii+1))/2,(y(jj)+y(jj+1))/2);
end
end
I=s*(b-a)*(d-c)/(m*n);
end

>> f = @ (x , y ) sin( x .* y )./ sqrt ( x + y );

>> a=0.5;b=1;c=0.5;d=2;
>> m=8;n=8;
>> I=Integralcenterpoint(f,a,b,c,d,n,m)
I=
0.3929
Note : Matlab has a built in command for double integrals on rectangles: integral2(f,a,b,c,d). Here is
an example:
>> f = @ (x , y ) sin( x .* y )./ sqrt ( x + y );

>> I = integral2 (f ,0.5 ,1 ,0.5 ,2)

I=

0.3924

The four corners method

20
The formula would be

Notice that the four-corners method coincides with applying the trapezoid rule in each direction. Thus it is
in fact a double trapezoid rule.
function wij = weightfourcorner(n,m)
if rem(m,2)~=0 || rem(n,2)~=0
error('not correct input/n&m even number ')
else
wij=2*ones(n+1,m+1);
wij(1,1)=1;
wij(n+1,1)=1;
wij(1,m+1)=1;
wij(n+1,m+1)=1;
wij(2:1:n,2:1:m)=4;
end
end

function I=fourcorner_I(f,a,b,c,d,n,m)
m;
n;
wij = weightfourcorner(n,m);
x=a:(b-a)/n:b;
y=c:(d-c)/m:d;
ss=0;
for i=1:1:(n+1)
for k=1:1:(m+1)
F=f(x(i),y(k)).*wij(i,k);
ss=ss+F;
end
end
I=((b-a)*(d-c)/(4*n*m)).*ss;
end

>> f = @ (x , y ) sin( x .* y )./ sqrt ( x + y );

>> a=0.5;b=1;c=0.5;d=2;
>> m=8;n=8;
>> wij = weightfourcorner(n,m)

21
 wij =

1 2 2 2 2 2 2 2 1
2 4 4 4 4 4 4 4 2
2 4 4 4 4 4 4 4 2
2 4 4 4 4 4 4 4 2
2 4 4 4 4 4 4 4 2
2 4 4 4 4 4 4 4 2
2 4 4 4 4 4 4 4 2
2 4 4 4 4 4 4 4 2
1 2 2 2 2 2 2 2 1

>> I=fourcorner_I(f,a,b,c,d,n,m)

I =
0.3913

The double Simpson method

Weight matrix of double Simpson method

The double Simpson formula

function w = mysimpweights(n)
if rem(n ,2) ~= 0
error('n must be even for Simpsons rule ')
end
w = 2* ones (n+1 ,1);
w(1) = 1; w(n+1) = 1;
w(2:2:n)=4;
end
function W = mydblsimpweights(m ,n)
if rem (m ,2)~=0 || rem (n ,2)~=0
error ('m and n must be even for Simpsons rule ')
end
u = mysimpweights (m);
v = mysimpweights (n);

22
W = u*v';
end
function I=doublesimpson_I(f,a,b,c,d,n,m)
m;
n;
wij = mydblsimpweights(m ,n);
x=a:(b-a)/n:b;
y=c:(d-c)/m:d;
ss=0;
for i=1:1:(n+1)
for k=1:1:(m+1)
F=f(x(i),y(k)).*wij(i,k);
ss=ss+F;
end
end
I=((b-a)*(d-c)/(9*n*m)).*ss;
end
>> f = @ (x , y ) sin( x .* y )./ sqrt ( x + y );
>> a=0.5;b=1;c=0.5;d=2;
>> m=8;n=8;
>> I=doublesimpson_I(f,a,b,c,d,n,m)

I=
0.3924

23
SECTION 5: Double Integrals for Non-rectangles
Example: find
When 1) T is the triangle with corners at (0, 0), (1, 0) and (0, 2).
2) T is the shaded region

Solution:
1) First case
>> f =@(x,y)sin(x.*y).^3.*(2*x+y <=2);
>> I = integral2 (f,0,1,0,2)

I=
0.0134
Or
1 𝑦=2−2𝑥
𝐼 = ∫0 ∫0 𝑠𝑖𝑛3 (𝑥𝑦)𝑑𝑦𝑑𝑥

>> f =@(x,y)sin(x.*y).^3;
>> xmin=0; xmax=1;
>> ymin=0; ymax=@(x)2-2*x;
>> I = integral2(f,xmin,xmax,ymin,ymax)
I=
0.0134

2) second case
>> f =@(x,y)sin(x.*y).^3.*(y<=3-x.^2 & y>=2*x);
>> I = integral2 (f,-3,1,-6,3)

I=

0.2059
Or
>> f =@(x,y)sin(x.*y).^3;
>> xmin=-3; xmax=1;
>> ymin=@(x)2*x; ymax=@(x)3-x.^2;
>> I = integral2(f,xmin,xmax,ymin,ymax)

I=
0.2059

24
Example for triple integral
1 1 𝜋 1 √1−𝑥 2 √1−𝑥 2 −𝑦 2
Find ∫−1 ∫0 ∫0 𝑦𝑠𝑖𝑛(𝑥) + 𝑧𝑐𝑜𝑠(𝑥)𝑑𝑧𝑑𝑦𝑑𝑥 and ∫−1 ∫−√1−𝑥2 ∫−√1−𝑥2 −𝑦2 𝑦𝑠𝑖𝑛(𝑥) + 𝑧𝑐𝑜𝑠(𝑥)𝑑𝑧𝑑𝑦𝑑𝑥

clear all
clc
fun = @(x,y,z) y.*sin(x)+z.*cos(x);
q = integral3(fun,0,pi,0,1,-1,1)
>> q =

2.0000
clear all
clc
fun = @(x,y,z) x.*cos(y) + x.^2.*cos(z)
xmin = -1;
xmax = 1;
ymin = @(x)-sqrt(1 - x.^2);
ymax = @(x) sqrt(1 - x.^2);
zmin = @(x,y)-sqrt(1 - x.^2 - y.^2);
zmax = @(x,y) sqrt(1 - x.^2 - y.^2);
q = integral3(fun,xmin,xmax,ymin,ymax,zmin,zmax)
q=

0.7796

25
SECTION 6: Numerical Differentiation
Approximating derivatives from data
Data from an experiment or a simulation:

Method of Numerical Formula Example

Differentiation for f(x)
𝒅𝒚 𝒙𝒆𝒙
Find 𝒅𝒙 |𝒙=𝟏 for 𝒚 = 𝒔𝒊𝒏(𝒙)
>> diff(x*exp(x)/sin(x),x)
>> subs(ans,1)
>>
double(ans)
ans =
4.3866

Forward Difference >> x=[1,1.1];

>> y=x.*exp(x)./sin(x);
>> yp=(y(2)-y(1))/(x(2)-x(1))
yp =
4.7759
Backward Difference >> x=[0.9,1];
>> y=x.*exp(x)./sin(x);
>> yp=(y(2)-y(1))/(x(2)-x(1))
yp =
4.0444
Central Difference >> x=[0.9,1,1.1];
>> y=x.*exp(x)./sin(x);
>> yp=(y(3)-y(1))/(x(3)-x(1))
yp =
4.4102

26
Errors of approximation

Partial Derivatives

27
Example:
clear all
clc
close all
x= 0:.1:5;
y=pi:.02*pi:2*pi;
f = @(x,y) y.*sin(x.*y)
for i=1:1:size(x,2)
for j=1:1:size(y,2)
fxy(i,j)=f(x(i),y(j));
end
end
fxy;
h=0.1;
k=0.02*pi;

28
for i=1:size(x,2)
for j=1:size(y,2)
if i>1&&i<size(x,2)&&j>1&&j<size(y,2)
dfxy(i,j)= (1/(4*h*k))*(fxy(i+1,j+1)-fxy(i+1,j-1)-fxy(i-1,j+1)+fxy(i-
1,j-1));

else
dfxy(i,j)=0;
end
end
end
dfxy=dfxy(2:size(x,2)-1,2:size(y,2)-1)
[X,Y] =meshgrid(0.1:.1:4.9
,1.02*pi:.02*pi:1.98*pi);
Z = dfxy
figure (1)
surf(X,Y,Z)%mesh(X,Y,Z)
%--------------------------------------
syms x y
fxy1 = y*sin(x*y)
ddfxy=diff(fxy1,x,y)
figure (2)
ezsurf(ddfxy,[0.1,4.9],[1.02*pi,1.98*pi])

29
The Main Sources of Error
Truncation Error
Truncation error is defined as the error caused directly by an approximation method.

Roundoff Error
In IEEE standard double precision (used by Matlab and most serious software), machine epsilon is 2−52 or
about 2.2×10−16.

>> (2^52+1) - 2^52

ans =
1
>> (2^53+1) - 2^53
ans =
0
Loss of Precision (also called Loss of Significance)
Bad Conditioning
Sources of error:
Truncation – the method is an approximation.
Roundoff – double precision arithmetic uses ≈ 15 significant digits.
Loss of precision – an amplification of roundoff error due to cancellations.
Bad conditioning – the problem is sensitive to input errors.
Error can build up after multiple calculations.

Example:
>> format long
>> format compact
>> f = @ ( x ) x ^2;
>> for i = 1:30;
h = 10^( - i )
df = ( f (1+ h ) - f (1))/ h;
relerr = (2 - df )/2
end
At first the relative error decreases since truncation error is reduced. Then loss of precision takes over
and the relative error increases to 1.

30
Example

31
SECTION 7: Differential Equations (ODEs /IVP)

32
Using Matlab to solve a system of ODE’s

In Matlab there are several commands that can be used to solve an initial value problem for a system of
differential equations. Each of these correspond to different solving methods. The standard one is
ode45, which uses the algorithm “Runge-Kutta 45”.

Example:
Solve 𝑦⃛ + 𝑒 𝑡 𝑦̈ + 𝑦̇ − 𝑦 = 𝑡𝑒 −𝑡 where 𝑡 ∈ [0,5], 𝑦(0) = 1, 𝑦̇ (0) = 0, 𝑦̈ (0) = −1
clear all
clc
dy = @(t,y)[y(2); y(3); y(1)-y(2)-exp(t)*y(3)+t*exp(-t)]
[T Y]= ode45(dy,[0 5] ,[1; 0;-1]);
plot(T, Y(:,1),'b')

33
Euler Methods
Numerical Solution of an IVP

34
The Euler Method

Solve 𝑦⃛ + 𝑒 𝑡 𝑦̈ + 𝑦̇ − 𝑦 = 𝑡𝑒 −𝑡 where 𝑡 ∈ [0,5], 𝑦(0) = 1, 𝑦̇ (0) = 0, 𝑦̈ (0) = −1

function[T,Y]= myeuler(f,tspan,y0,n) dy = @(t,y)[y(2); y(3); y(1)-y(2)-
a = tspan (1); b = tspan (2); exp(t)*y(3)+t*exp(-t)];
h = (b-a)/n; >> [T Y] = myeuler (dy ,[0 5] ,[1; 0;-1] ,200);
t = a; T = a; >> theta = Y (: ,1);
y = y0; >> plot(T , theta )
Y = y0';
for i = 1:n
y = y + h*f(t,y);
t = a + i*h;
T = [T; t];
Y = [Y; y'];
end
end

35
The Modified Euler Method

function [T , Y] = >> dy = @(t,y)[y(2); y(3); y(1)-y(2)-exp(t)*y(3)+t*exp(-t)];

mymodeuler(f,tspan,y0,n) >> [T Y] = mymodeuler (dy ,[0 5] ,[1; 0;-1] ,200);
a = tspan(1); b = >> theta = Y (: ,1);
tspan(2); >> plot(T , theta )
h = (b-a)/n;
t = a; T = a;
y = y0;
Y = y0';
for i = 1:n
k1 = h*f(t,y);
k2 = h*f(t+h,y+k1);
y = y + .5*(k1+k2);
t = a + i*h;
T = [T; t];
Y = [Y; y'];
end
end

36
Fourth Order Runge-Kutta

The most famous of all IVP methods is the classic Runge-Kutta method of order 4:

Thus it is not easy beforehand to tell how small h should be to get the error within a given tolerance. For instance, if
the true solution oscillates very rapidly, we will obviously need a smaller step size than for a solution that is nearly
constant. How can a program then choose h small enough to produce the required accuracy? We also do not wish to
make h much smaller than necessary, since that would increase the number of steps. To accomplish this a program
tries an h and tests to see if that h is small enough. If not it tries again with a smaller h. If it is too small, it accepts that
step, but on the next step it tries a larger h. This process is called variable step size.

The Runge-Kutta 45 method, which is used in ode45, is an embedded method. In the RK45, the function f is evaluated
at 5 different points. These are used to make a 5th order estimate yi+1. At the same time, 4 of the 5 values are used to
also get a 4th order estimate. If the 4th order and 5th order estimates are close, then we can conclude that they are
accurate. If there is a large discrepency, then we can conclude that they are not accurate and a smaller h should be
used.

To see variable step size in action, we will define and solve two different ODEs and solve them on the same interval.
Create this script and run it:

37
% illustrates variable step size in RK45
dy1 = @(t , y) [-y(2); y(1)]; % create two ODE IVPs
dy2 = @(t , y) [-5*y(2);5*y(1)];
[T1 Y1]=ode45(dy1,[0 20],[1;0]); % solve with ode45
[T2 Y2]=ode45(dy2,[0 20],[1;0]);
y1 = Y1 (: ,1); % extract position variables
y2 = Y2 (: ,1);
plot ( T1 , y1 , 'b') % plot both together
hold on
plot ( T2 , y2 , 'r')
size ( T1 ) % print number of steps used
size ( T2 )
hold off

38
clear all
clc
%problem with the Euler method
%x''+x=0 y1=x y2=x' y1'=y2 y2'=-y1
% solution y1=cos(t) y2=sin(t)
xx=3; %xx=1 euler method, xx=2 modefied euler method
if xx==1
dy = @(t,y)[y(2); -y(1)]
[T Y] = myeuler (dy ,[0 2*pi] ,[1;0] ,500)
y1 = Y(: ,1);
y2 = Y(: ,2);
plot(y1 ,y2,'b')
hold on
else if xx==2
%The Modified Euler Method
dy = @(t,y)[y(2); -y(1)]
[T Ym] = mymodeuler (dy ,[0 2*pi] ,[1;0] ,20)
ym1 = Ym(: ,1);
ym2 = Ym(: ,2);
plot(ym1 ,ym2,'r')
else
dy = @(t,y)[y(2); -y(1)]
[T YR] =ode45(dy ,[0 2*pi] ,[1;0]);
yR1 = YR(: ,1);
yR2 = YR(: ,2);
plot(yR1 ,yR2,'g')
end
end

39
SECTION 8: Differential Equations (ODE/BVP and and Finite Differences)
Application Examples
Steady State Heat and Diffusion

Beam With Tension

Consider a simply supported beam with modulus of elasticity E, moment of inertia I, a uniform load w, and end tension
T. If y(x) denotes the deflection at each point x in the beam, then y(x) satisfies the differential equation

40
Finite Difference Method – Linear ODE
clear all
clc
%Finite Dierence Method { Linear ODE BVP}
%y''+y'-y=x y(0)=y(1)=0;
t=0:0.01:1;
y=exp(-t/2).*(cosh((5^(1/2)*t)/2) +
(13133*5^(1/2)*sinh((5^(1/2)*t)/2))/25000) - t - 1;%exact
plot(t,y,'b')
hold on
%({y(i+1)-2y(i)+y(i-1)}/h^2)+(y(i+1)-y(i-1))/(2h) -y(i)=x(i)
%x= 0 0.2 0.4 0.6 0.8 1
%i= 0 1 2 3 4 5
%y= 0 y1 y2 y3 y4 0
%27.5y2-51y1=0.2 -------1)
%27.5y3-51y2+22.5y1=0.4 ------2)
%27.5y4-51y3+22.5y2=0.6 -------3)
%-51y4+22.5y3=0.8 --------4)
A=[-51 27.5 0 0;22.5 -51 27.5 0;0 22.5 -51 27.5;0 0 22.5 -51];
B=[0.2 0.4 0.6 0.8]';
y1=A\B;
x = 0:0.2:1;
yx = [0 ; y1 ; 0]
plot(x,yx,'ro-')

Exercise:

41
Finite Difference Method – Nonlinear ODE
If we again consider the heat in a metal bar of length L, but this time consider the effect of radiation as
well as conduction, then the steady state equation has the form

where ub is the temperature of the background, d incorporates a coefficient of radiation and g(x) is the
heat source.
If we again replace the continuous problem by its discrete approximation then we get

This equation is nonlinear in the unknowns, thus we no longer have a system of linear equations to solve,
but a system of nonlinear equations. One way to solve these equations would be by the multivariable
Newton method.
Example: Solve the a bove finite difference equations with the boundary conditions u(0) = 0 and u(L) = 0.
We let L = 4, n = 4, d = 1, and g(x) = sin(πx/4).
clear all
clc
syms u1 u2 u3
f1=u2-2*u1-(u1^4-0.5^4)+sin(pi/4);
f2=u3-2*u2+u1-(u2^4-0.5^4)+sin(2*pi/4);
f3=-2*u3+u2-(u3^4-0.5^4)+sin(3*pi/4);
f=[f1;f2;f3];
Df=[diff(f1,u1) diff(f1,u2) diff(f1,u3)
diff(f2,u1) diff(f2,u2) diff(f2,u3)
diff(f3,u1) diff(f3,u2) diff(f3,u3)];
n=20;
u123 = [0.75;0.9;0.75];
for i= 1:n
Du123=-double(subs(Df,[u1,u2,u3],u123')\subs(f,[u1,u2,u3],u123'));
u123 = u123 + Du123;
f=double(subs(f,[u1;u2;u3],u123));
end
u123;
f;
x = 0:4;
u = [0 ; u123 ; 0];
plot(x,u,'ro-')
hold on

42
Relaxation Method for Nonlinear Finite Differences

43
Example: Solve the a bove finite difference equations with the boundary conditions u(0) = 0 and u(L) = 0.
We let L = 4, n = 4, d = 1, and g(x) = sin(πx/4).
L = 4;
n = 4;
h = L/n ;
hh = (h^2)/3;
u0 = 0;
uL = 0;
ub = .5;
ub4 = ub ^4;
x = 0: h : L ;
g = sin(pi* x/4);
u = zeros (1,n +1);
steps = 4;
u (1)= u0 ;
u (n+1)= uL ;
for j = 1: steps
u (2:n)=(u(3:n+1)+u(2:n)+u(1:n-1))/3 + hh *(-u(2:n).^4+ ub4 + g(2:n));
end
plot(x,u)

If you run this program with the given n and steps the result will not seem reasonable. We can plot the
initial guess by adding the command plot(x,u) right before the for loop. We can also plot successive
iterations by moving the last plot(x,u) before the end. Now we can experiment and see if the iteration is
converging. Try various values of steps and n to produce a good plot. You will notice that this method
converges quite slowly. In particular, as we increase n, we need to increase steps like n 2 , i.e. if n is large
then steps needs to be really large.

44
SECTION 9: Parabolic PDEs

Introduction of Eigenvalues and Eigenvectors

Suppose that A is a square (n × n) matrix. We say that a nonzero vector v is an eigenvector and a number
λ is its eigenvalue if
Av = λv → det(A − λI) = 0
Geometrically this means that Av is in the same direction as v.
>> A = pascal (4)
A=
1 1 1 1
1 2 3 4
1 3 6 10
1 4 10 20
>> [ v e ] = eig( A )→ function for computing Eigenvalues and Eigenvectors
v=
0.3087 -0.7873 0.5304 0.0602
-0.7231 0.1632 0.6403 0.2012
0.5946 0.5321 0.3918 0.4581
-0.1684 -0.2654 -0.3939 0.8638
e=

0.0380 0 0 0
0 0.4538 0 0
0 0 2.2034 0
0 0 0 26.3047

>> v1 = v (: ,1) Note: A symmetric matrix is called positive definite if for

v1 = all vectors v ≠ 0 the following holds: Av · v > 0.
0.3087 Geometrically, A does not rotate any vector by more than
-0.7231 π/2. In summary:
0.5946 • If A is symmetric then its eigenvalues are real.
-0.1684 • If A is symmetric positive definite, then its eigenvalues
are positive numbers
>> A * v1
ans =
0.0117
-0.0275
0.0226
-0.0064

>> e (1 ,1)* v1
ans =

0.0117
-0.0275
0.0226
-0.0064

45
Explicit Method
Heat Flow and Diffusion:

In three dimensions, the equation that governs both of these processes is the heat/diffusion equation
,
where c is the coefficient of conduction or diffusion, and ∆u(x, y, z) = uxx +uyy +uzz. The symbol ∆ in this
context is called the Laplacian.

In some problems the z dimension is irrelevent, either because the object in question is very thin, or u does
not change in the z direction. In this case the equation

Explicit Method Finite Differences

Here we have used the forward difference for ut and the central difference for uxx. This equation can be
solved for ui,j+1 to produce

46
Initial Condition

The value at grid point (i, j+1) depends on its previous value and the previous values of its nearest
neighbors.

Boundary Conditions
at x = 0 and x = L. One possibility is fixed boundary conditions, which we can implement just as we did
for the ODE boundary value problem.

In a heat problem, g1 and g2 would represent heating or cooling applied to the ends. These are easily
implemented in a program by letting u0,j = g1(tj ) and um,j = g2(tj ).

% solve ut = cuxx for 0≤x≤L, 0≤t≤T

% BC: u(0,t) = g1(t); u(L,t)=g2(t)
% IC: u(x,0) = f(x)
% Inputs:
% f -- function for IC
% g1,g2 -- functions for BC
% L -- length of rod
% T -- length of time interval
% m -- number of subintervals for x
% n -- number of subintervals for t
% c -- rate constant in equation
% Outputs:
% t -- vector of time points
% x -- vector of x points
% u -- matrix of the solution, u(i,j)~=u(x(i),t(j))
function [t x u] = myheat(f,g1,g2,L,T,m,n,c)
h = L/m; k = T/n; % set space and time step sizes
r = c*k/h^2; rr = 1-2*r;
x = linspace(0,L,m+1); % set space discretization
t = linspace(0,T,n+1); % set time discretization

u = zeros(m+1,n+1); %Set up the matrix for u:

u(:,1) = f(x); % evaluate initial conditions
u(1,:)= g1(t); u(m+1,:)= g2(t); % evaluate boundary conditions
for j = 1:n

47
 u(2:m,j+1) = r*u(1:m-1,j) + rr*u(2:m,j) + r*u(3:m+1,j)
end
%figure (1)
%mesh(x,t,u')
figure (2)
surf(x,t,u')
Example 1:
>> m=10; n=10; c=1;L = 20; T = 20; f=@(x) 0.5*x; g1 =@(t)0; g2= @(t) cos(t);
>> [t x u] = myheat(f,g1,g2,L,T,m,n,c);

Example 2:
>> m=20;n=100;c=0.25;L=2*pi;T=20;g1=@(t)sin(t);g2=@(t)0;f=@(x)-sin(x/4);
>> [t x u] = myheat(f,g1,g2,L,T,m,n,c);

Example 3:
>> m=10;n=20;c=1;L=1;T=1/3;g1=@(t)0;g2=@(t)0;f=@(x)6*sin(pi*x);
>> myheat(f,g1,g2,L,T,m,n,c);
Exact solution: u= @(x,t) 6*sin(pi*x)*exp(-(pi^2)*t);
ezsurf(u,[0,1],[0,1/3]);

Numerical solution Exact solution

48
 Example 4:
>> m=25;n=35;c=0.25;L=pi;T=1;g1=@(t)0;g2=@(t)0;f=@(x)3*sin(5*x/2);
>> myheat(f,g1,g2,L,T,m,n,c);
Exact solution: u= @(x,t) 3*sin(5*x/2)*exp(-((5*0.25/2)^2)*t);
ezsurf(u,[0,pi],[0,1]);

Numerical solution Exact solution

Solution Instability for the Explicit Method

As we saw in experiments using myheat.m, the solution can become unbounded unless the time steps are
small. In this section we consider why.

Writing the Difference Equations in Matrix Form:

49
The problem with this is the term AkEj . This term is exactly what we would do in the power method for
finding the eigenvalue of A with the largest absolute value. If the matrix A has eigenvalues with absolute
value greater than 1, then this term will grow exponentially. The following Figure shows the largest absolute
value of an eigenvalue of A as a function of the parameter r for various sizes of the matrix A. As you can
see, for r > 1/2 the largest absolute eigenvalue grows rapidly for any m and quickly becomes greater than
1.

Maximum absolute eigenvalue as a function of r for the matrix A from the explicit method for the heat
equation calculated for matrices A of sizes m = 2 . . . 10. Whenever the maximum absolute eigenvalue is
greater than 1 the method is unstable, i.e. errors grow exponentially with each step. When using the explicit
method r < 1/2 is a safe choice.

Recall that r = ck/h2 . Since this must be less than 1/2, we have

50
Test this using m = 6 and r = .4, .6. Check the eigenvalues and eigenvectors of the resulting matrices:

>> A = myexpmatrix_A (6 ,.6)

A=

-0.2000 0.6000 0 0 0 0
0.6000 -0.2000 0.6000 0 0 0
0 0.6000 -0.2000 0.6000 0 0
0 0 0.6000 -0.2000 0.6000 0
0 0 0 0.6000 -0.2000 0.6000
0 0 0 0 0.6000 -0.2000

>> [v e] =eig(A)

v=

0.2319 -0.4179 -0.5211 0.5211 -0.4179 0.2319

-0.4179 0.5211 0.2319 0.2319 -0.5211 0.4179
0.5211 -0.2319 0.4179 -0.4179 -0.2319 0.5211
-0.5211 -0.2319 -0.4179 -0.4179 0.2319 0.5211
0.4179 0.5211 -0.2319 0.2319 0.5211 0.4179
-0.2319 -0.4179 0.5211 0.5211 0.4179 0.2319

e=

-1.2812 0 0 0 0 0
0 -0.9482 0 0 0 0
0 0 -0.4670 0 0 0
0 0 0 0.0670 0 0
0 0 0 0 0.5482 0
0 0 0 0 0 0.8812

Cheke
Ak=vekv-1
>> A = myexpmatrix_A (6 ,.4)

A=

51
 0.2000 0.4000 0 0 0 0
0.4000 0.2000 0.4000 0 0 0
0 0.4000 0.2000 0.4000 0 0
0 0 0.4000 0.2000 0.4000 0
0 0 0 0.4000 0.2000 0.4000
0 0 0 0 0.4000 0.2000

>> [v e] =eig(A)

v=

0.2319 -0.4179 -0.5211 0.5211 -0.4179 0.2319

-0.4179 0.5211 0.2319 0.2319 -0.5211 0.4179
0.5211 -0.2319 0.4179 -0.4179 -0.2319 0.5211
-0.5211 -0.2319 -0.4179 -0.4179 0.2319 0.5211
0.4179 0.5211 -0.2319 0.2319 0.5211 0.4179
-0.2319 -0.4179 0.5211 0.5211 0.4179 0.2319

e=

-0.5208 0 0 0 0 0
0 -0.2988 0 0 0 0
0 0 0.0220 0 0 0
0 0 0 0.3780 0 0
0 0 0 0 0.6988 0
0 0 0 0 0 0.9208

52
The Implicit Difference Equations

53
Section 10: Finite Difference Method for Elliptic PDEs

where:

The boundary conditions are introduced by:

54
Matlab program for solving poisson equation:
clear all;close all;clc;
% mypoisson
% Poisson on a Rectangle
% Solves u_xx + u_yy = f(x,y)
% with boundary conditions:
% u(a,y) = u(x,c) = 1
% u(b,y) = 1 + 5/8 y
% u(x,d) = 1 + 1/2 x
% with f(x,y) = .5 sin(x)sin(y)
% Define the rectangle and nodes
a = 0;b = 10;c = 0;d = 8;
m = 15;n = 10;
h = (b-a)/m;
k = (d-c)/n;
x = a:h:b;
y = c:k:d;
% Assign number of iterations and calculate constants
maxit = 0;
c1 = h^2/(2*(h^2 + k^2));
c2 = k^2/(2*(h^2 + k^2));
c3 = k^2*c1;
% Assign initial values and boundary conditions
u = zeros(m+1,n+1);
u(1,:) = 1;
u(:,1) = 1;
u(:,n+1) = 1 + 5*sin(x(:)*pi/b/2);
u(m+1,:) = 1 + (5/8)*y(:);
% Assign values of f(x,y) in the interior
f = zeros(m+1,n+1);
for i = 2:m
f(i,:) = .5*sin(x(i)).*sin(y(:));
end
% Iterate values of u in the interior
for j = 1:maxit
u(2:m,2:n) = c1*(u(2:m,3:n+1)+u(2:m,1:n-1)) + c2*(u(3:m+1,2:n)+u(1:m-
1,2:n)) - c3*f(2:m,2:n);
end
% Plot the result
figure (1)
mesh(x,y,u')
figure (2)
surf(x,y,u')

55
You will notice that maxit is set to 0. Thus, the program will not do any iteration but will plot the initial
guess. The initial guess in this case consists of the proper boundary values at the edges and zero everywhere
in the interior. To see the solution evolve, gradually increase maxit.

maxit = 0

maxit = 1

maxit = 2

maxit = 5

maxit = 100

56