Dear Student, Here you can find MATLAB code of some numerical techniques. You can modify the code if you find error in execution!!
BISECTION Algorithm:
>> x = [-3:0.1:3];
>> y = fnroot(x);
>> plot(x,y);
>> grid on;
>> ylabel('x^3 - x -1');
>> xlabel('x');
>> title('Graph of f(x) cuts x-axis');
>> ylabel('f(x)=x^3 - x -1');
>> x = [-3:0.1:3];
>> y = fnroot(x);
>> plot(x,y);
>> grid on;
>> ylabel('x^3 - x -1');
>> xlabel('x');
>> title('Graph of f(x) cuts x-axis');
>> ylabel('f(x)=x^3 - x -1');
function y = fnroot(x)
y = x.^3-x-1;
end
function [ c,err ] = bisectf(f,a,b,E)
fa = feval(f,a);
fb = feval(f,b);
for i=1 : 30
c(i) =(a+b)/2;
fc = feval(f,c(i));
err(i)= abs(a-b);
if abs(fc) < E
break
else if fa * fc < 0
b = c(i);
fb = fc;
else
a = c(i);
fa = fc;
end
end
end
end
>> [c ,err] = bisectf('fnroot',0,2,0.005);
>> [c' err']
ans =
1.0000 2.0000
1.5000 1.0000
1.2500 0.5000
1.3750 0.2500
1.3125 0.1250
1.3438 0.0625
1.3281 0.0313
1.3203 0.0156
1.3242 0.0078
NEWTON RAPHSON Algorithm
y = x.^3-x-1;
end
function y = df(x)
y = 3*x*x-1;
end
function [ n,c ] = ntr(f,df,a,E)
fa = feval(f,a);
dfa = feval(df,a);
for i=1 : 30
c(i) = a - fa/dfa;
fa = feval(f,c(i));
n(i) = i;
if abs(fa) < E
break
else if abs(a - c(i))< E
break
else
dfa = feval(df,c(i));
a = c(i);
end
end
end
end
>> [n,root]=ntr('fnroot','df',1,0.0005);
>> [n' root']
ans =
1.0000 1.5000
2.0000 1.3478
3.0000 1.3252
4.0000 1.3247
NEWTON'S FORWARD DIFFERENCE Interpolation Algorithm (File Name Forward_D.m)
x = input('Equidistant values of x:');
y = input('Values of y for each x:');
n = length(x);
d = zeros(n);
for i=1:n
d(i,1)=y(i);
end
for j=2:n
for i=1:n-(j-1)
d(i,j)=d(i+1,j-1)-d(i,j-1);
end
end
D = [x',d];
fprintf('The Forward Difference Table\n');
disp(D);
%disp(d);
a = input('Enter value of x to find y:');
h = x(2)-x(1);
product =1;
sum = d(1,1);
p = (a-x(1))/h;
for j=2:n
product = product * (p-(j-2))/(j-1);
sum = sum + product*d(1,j);
end
>> Forward_D
Equidistant values of x:[0 1 2 3]
Values of y for each x:[1 0 1 10]
The Forward Difference Table
0 1 -1 2 6
1 0 1 8 0
2 1 9 0 0
3 10 0 0 0
Enter value of x to find y:4
The Interpolating value y(4) is 33.000000
Solution to Linear System of Equation AX = b
Gauss Jacobi Algorithm
>> A = [10 -2 -1 -1;-2 10 -1 -1;-1 -1 10 -2;-1 -1 -2 10];
>> b = [3;15;27;-9];
>> A
A =
10 -2 -1 -1
-2 10 -1 -1
-1 -1 10 -2
-1 -1 -2 10
>> b
b =
3
15
27
-9
>>
function x = jacobinew( A,b,M,E)
n = length(b);
fprintf('The Solution to sytem of equation AX = b is \n');
for i=1:n
if i==1
B(i,:)=[A(i,2),A(i,3:n)]*(-1/A(i,i));
else if i==n
B(i,:)=[A(i,1:n-1)]*(-1/A(i,i));
else
B(i,:)=[A(i,1:i-1),A(i,i+1:n)]*(-1/A(i,i));
end
end
end
for i=1:n
C(i)=b(i)*(1/A(i,i));
end
C = C';
x = zeros(n,1);
u = zeros(n,1);
y = zeros(n);
temp = 0;
for j=2:M
for i=1:n
if i==1
y(i,j)= B(i,:)*u(2:n)+C(i);
x(i)=B(i,:)*u(2:n)+C(i);
else if i==n
y(i,j)=B(i,:)*u(1:n-1)+C(i);
x(i)=B(i,:)*u(1:n-1)+C(i);
else
w = [u(1:i-1);u(i+1:n)];
y(i,j)=B(i,:)*w+C(i);
x(i)=B(i,:)*w+C(i);
end
end
end
u = x;
disp(x')
for i=1:n
if abs(y(i,j)-y(i,j-1)) > E
temp = temp + 1;
end
end
if temp ==0
break;
end
temp = 0;
end
end
The Solution to sytem of equation AX = b is
0.3000 1.5000 2.7000 -0.9000
0.7800 1.7400 2.7000 -0.1800
0.9000 1.9080 2.9160 -0.1080
0.9624 1.9608 2.9592 -0.0360
0.9845 1.9848 2.9851 -0.0158
0.9939 1.9938 2.9938 -0.0060
0.9975 1.9975 2.9976 -0.0025
0.9990 1.9990 2.9990 -0.0010
0.9996 1.9996 2.9996 -0.0004
>> x = jacobinew(A,b,15,0.0001);
The Solution to sytem of equation AX = b is
0.3000 1.5000 2.7000 -0.9000
0.7800 1.7400 2.7000 -0.1800
0.9000 1.9080 2.9160 -0.1080
0.9624 1.9608 2.9592 -0.0360
0.9845 1.9848 2.9851 -0.0158
0.9939 1.9938 2.9938 -0.0060
0.9975 1.9975 2.9976 -0.0025
0.9990 1.9990 2.9990 -0.0010
0.9996 1.9996 2.9996 -0.0004
0.9998 1.9998 2.9998 -0.0002
0.9999 1.9999 2.9999 -0.0001
>> x = jacobinew(A,b,15,0.00001);
The Solution to sytem of equation AX = b is
0.3000 1.5000 2.7000 -0.9000
0.7800 1.7400 2.7000 -0.1800
0.9000 1.9080 2.9160 -0.1080
0.9624 1.9608 2.9592 -0.0360
0.9845 1.9848 2.9851 -0.0158
0.9939 1.9938 2.9938 -0.0060
0.9975 1.9975 2.9976 -0.0025
0.9990 1.9990 2.9990 -0.0010
0.9996 1.9996 2.9996 -0.0004
0.9998 1.9998 2.9998 -0.0002
0.9999 1.9999 2.9999 -0.0001
1.0000 2.0000 3.0000 -0.0000
1.0000 2.0000 3.0000 -0.0000
1.0000 2.0000 3.0000 -0.0000
>>
Gauss Seidel Algorithm
function x = gseidel(A,b,M,E)
n = length(b);
fprintf('The Solution to system of equation AX = b is \n');
for i=1:n
if i==1
B(i,:)=[A(i,2:n)]*(-1/A(i,i));
else if i==n
B(i,:)=[A(i,1:n-1)]*(-1/A(i,i));
else
B(i,:)=[A(i,1:i-1),A(i,i+1:n)]*(-1/A(i,i));
end
end
end
for i=1:n
C(i)=b(i)/A(i,i);
end
C = C';
x =zeros(n,1);
y = zeros(n);
temp = 0;
for j= 2:M
for i=1:n
if i==1
x(i) = B(i,:)*x(2:n)+C(i);
else if i==n
x(i)=B(i,:)*x(1:n-1)+C(i);
else
x(i)=B(i,:)*[x(1:i-1);x(i+1:n)]+C(i);
end
end
y(i,j) = x(i);
end
for i=1:n
if abs(y(i,j)-y(i,j-1)) > E
temp = temp + 1;
end
end
if temp ==0
break;
end
disp(x');
temp = 0;
end
end
>> x = gseidel(A,b,15,0.001);
The Solution to system of equation AX = b is
0.3000 1.5600 2.8860 -0.1368
0.8869 1.9523 2.9566 -0.0248
0.9836 1.9899 2.9924 -0.0042
0.9968 1.9982 2.9987 -0.0008
0.9994 1.9997 2.9998 -0.0001
>> x = gseidel(A,b,15,0.00001);
The Solution to system of equation AX = b is
0.3000 1.5600 2.8860 -0.1368
0.8869 1.9523 2.9566 -0.0248
0.9836 1.9899 2.9924 -0.0042
0.9968 1.9982 2.9987 -0.0008
0.9994 1.9997 2.9998 -0.0001
0.9999 1.9999 3.0000 -0.0000
1.0000 2.0000 3.0000 -0.0000
1.0000 2.0000 3.0000 -0.0000
>>
Numerical Integration Algorithm
function t = f(x)
t = (1+x).^(-1);
end
function I = simpson13(y,n,h)
sum13 = y(1)+y(n+1);
for i=2:n
if rem(i,2)==0
sum13 = sum13 + 4*y(i);
else
sum13 = sum13 + 2*y(i);
end
I = (h/3)*sum13;
end
function I = Simpson38(y,n,h)
sum38 = y(1)+y(n+1);
for i=2:n
if rem(i-1,3)==0
sum38 = sum38 + 2*y(i);
else
sum38 = sum38 + 3*y(i);
end
end
I = (3*h/8)*sum38;
end
function I = trapz(y,n,h)
sumtr = y(1)+y(n+1);
for i=2:n
sumtr = sumtr + 2*y(i);
end
I = (h/2)*sumtr;
end
function Itg = NumInteg(f)
a = input('Enter lower limit of integration:');
b = input('Enter Upper limit of integration:');
n = input('Number of Intervals:');
h = (b-a)/n;
x = a:h:b;
y = feval(f,x);
fprintf('Table values of (x,y) is\n');
[x' y']
fprintf('The Function is f(x) = 1/(1+x)\n');
if mod(n,2)== 0
fprintf('Simpson 1/3 Rule:\n');
Itg = simpson13(y,n,h);
elseif mod(n,3) == 0
fprintf('Simpson 3/8 Rule:\n');
Itg = Simpson38(y,n,h);
else
fprintf('Trapezoidal Rule:\n');
Itg = trapz(y,n,h);
end
>> I = NumInteg('f')
Enter lower limit of integration:0
Enter Upper limit of integration:1
Number of Intervals:4
Table values of (x,y) is
ans =
0 1.0000
0.2500 0.8000
0.5000 0.6667
0.7500 0.5714
1.0000 0.5000
The Function is f(x) = 1/(1+x)
Simpson 1/3 Rule:
I =
0.6933
>> I = NumInteg('f')
Enter lower limit of integration:0
Enter Upper limit of integration:1
Number of Intervals:6
Table values of (x,y) is
ans =
0 1.0000
0.1667 0.8571
0.3333 0.7500
0.5000 0.6667
0.6667 0.6000
0.8333 0.5455
1.0000 0.5000
The Function is f(x) = 1/(1+x)
Simpson 1/3 Rule:
I =
0.6932
>> I = NumInteg('f')
Enter lower limit of integration:0
Enter Upper limit of integration:1
Number of Intervals:7
Table values of (x,y) is
ans =
0 1.0000
0.1429 0.8750
0.2857 0.7778
0.4286 0.7000
0.5714 0.6364
0.7143 0.5833
0.8571 0.5385
1.0000 0.5000
The Function is f(x) = 1/(1+x)
Trapezoidal Rule:
I =
0.6944
>>
Dear Student, you can download here Simple Algorithms of some numerical techniques that you know. You can use them to write a simple C/C++ program or MATLAB Program.
No comments:
Post a Comment