//Simpson's (3/8)th Rule
deff('y=f(x)','y=x^4')
a=input("Enter Lower Limit: ")
b=input("Enter Upper Limit: ")
n=input("Enter number of sum intervals: ")
h=(b-a)/n
add1=0
add2=0
add3=0
for i=0:n
x=a+i*h
y=f(x)
disp([x y])
if (i==0)|(i==n) then
add1=add1+y
else if (modulo(i,3)==0) then
add2=add2+y
else
add3=add3+y
end
end
end
I=((3*h)/8)*(add1+2*add2+3*add3)
disp(I,"Integration by Simpsons (3/8)th Rule is:")
Output
-->exec('D:\Scilab prog by me\Simpson's (3-8)th Rule 2.sce', -1)
Enter Lower Limit: -3
Enter Upper Limit: 3
Enter number of sum intervals: 12
- 3. 81.
- 2.5 39.0625
- 2. 16.
- 1.5 5.0625
- 1. 1.
- 0.5 0.0625
0. 0.
0.5 0.0625
1. 1.
1.5 5.0625
2. 16.
2.5 39.0625
3. 81.
Integration by Simpsons (3/8)th Rule is:
97.3125
deff('y=f(x)','y=x^4')
a=input("Enter Lower Limit: ")
b=input("Enter Upper Limit: ")
n=input("Enter number of sum intervals: ")
h=(b-a)/n
add1=0
add2=0
add3=0
for i=0:n
x=a+i*h
y=f(x)
disp([x y])
if (i==0)|(i==n) then
add1=add1+y
else if (modulo(i,3)==0) then
add2=add2+y
else
add3=add3+y
end
end
end
I=((3*h)/8)*(add1+2*add2+3*add3)
disp(I,"Integration by Simpsons (3/8)th Rule is:")
Output
-->exec('D:\Scilab prog by me\Simpson's (3-8)th Rule 2.sce', -1)
Enter Lower Limit: -3
Enter Upper Limit: 3
Enter number of sum intervals: 12
- 3. 81.
- 2.5 39.0625
- 2. 16.
- 1.5 5.0625
- 1. 1.
- 0.5 0.0625
0. 0.
0.5 0.0625
1. 1.
1.5 5.0625
2. 16.
2.5 39.0625
3. 81.
Integration by Simpsons (3/8)th Rule is:
97.3125
hey when i put 9 instead of 12 i m getting wrong answers
ReplyDeleteclc;clear;
ReplyDelete//using Trapezoidal rule, evaluate integral from 0 to 1 e^x , n=4
function[i]=trapezoidal(a,b,n,f)
h=(b-a)/n;
x=(a:h:b);
y=f(x);
m=length(y);
i=y(1)+y(m);
for j=2:m-1
i=i+2*y(j);
end
i=(h/2)*i;
disp('ans=',i)
return(i);
endfunction
deff('[y]=f(x)','y=exp(x)');
trapezoidal(0,1,4,f)
clc;clear;close;
ReplyDelete//using Simpson's 1/3rd rule, evaluate integral from 0 to 4 e^x , n=4
function[i]=simpsons13(a,b,n,f)
h=(b-a)/n;
x=(a:h:b);
y=f(x);
m=length(y);
i=y(1)+y(m);
for j=2:m-1
if (modulo(j,2)==0) then
i=i+2*y(j);
else
i=i+4*y(j);
end
end
i=(h/3)*i;
disp('ans=',i);
return(i);
endfunction
deff('[y]=f(x)','y=exp(x)');
simpsons13(0,4,4,f)
// Practical no. 5
ReplyDelete//Q 1 solve following DE to find y(0.2) using Euler's method.
//dy/dx = y-xy+x, y(0)=1.
function[y0]=eular(x0,y0,h,yest,f)
n=(yest-x0)/h;
for i=1:n
y0=y0+f(x0,y0)*h;
x0=x0+h;
disp(y0);
end
endfunction
deff('[y0]=f(a,b)','y=b-a*b+a');
eular(0,1,0.2,1,f)