Thread: Simpsons 3.8 rule, error help matlab

1. Simpsons 3.8 rule, error help matlab

hi this is my code

Code:
function [y]=BJ(n,x)
N=100;
b=pi;
a=0;
h=(b-a)/N;
% Remeber that y(1) is the first element, not y(0).
S=0;
for i=0:N-1
tj=i*h;
delta=a:b;
y(1)=cos((n*tj)-x*sin(tj));
y(2)=cos((n*(tj+h/3))-x*sin(tj+h/3));
y(3)=cos((n*(tj+(2*h/3)))-x*sin(tj+(2*h/3)));
y(4)=cos((n*(tj+h))-x*sin(tj+h));
y4(delta)=((x^4)*((cos(delta))^4)-4*n*(x^3)*((cos(delta))^3)+(4+6*(n^2))*(x^2)*((cos(delta))^2)-4*x*n*(1+(n^2))*cos(delta)+(n^4)-3*(x^2)*((sin(delta))^2))*cos(n*delta-x*sin(delta))+6*sin(n*delta- x*sin(delta))*((n^2)-2*x*cos(delta)*n+(x^2)*((cos(delta))^2)+(1/6))*x*sin(delta);
S = S+((h/8)*y(1))+((3*h/8)*y(2))+((3*h/8)*y(3))+((h/8)*y(4))-((3*(h^5)/80)*(y4(delta)));
end
y=1/pi*S;
end
can anyone fix it for me, to include the error... i got it right without it, but for some reason putting the error in, stuffs it up, and cant find where

it says something about the matrix must be square????

2. also does anyone know why, my 3/8 rule is different to the one on the net??? mine is h/8(f0+3f1+3f2+f3) and theirs is 3h/8(f0+3f1+3f2+f3) and my code works and is right according to matlab's inbuilt besselj function?

3. Originally Posted by fredrick08
hi this is my code

Code:
function [y]=BJ(n,x)
N=100;
b=pi;
a=0;
h=(b-a)/N;
% Remeber that y(1) is the first element, not y(0).
S=0;
for i=0:N-1
tj=i*h;
delta=a:b;
y(1)=cos((n*tj)-x*sin(tj));
y(2)=cos((n*(tj+h/3))-x*sin(tj+h/3));
y(3)=cos((n*(tj+(2*h/3)))-x*sin(tj+(2*h/3)));
y(4)=cos((n*(tj+h))-x*sin(tj+h));
y4(delta)=((x^4)*((cos(delta))^4)-4*n*(x^3)*((cos(delta))^3)+(4+6*(n^2))*(x^2)*((cos(delta))^2)-4*x*n*(1+(n^2))*cos(delta)+(n^4)-3*(x^2)*((sin(delta))^2))*cos(n*delta-x*sin(delta))+6*sin(n*delta- x*sin(delta))*((n^2)-2*x*cos(delta)*n+(x^2)*((cos(delta))^2)+(1/6))*x*sin(delta);
S = S+((h/8)*y(1))+((3*h/8)*y(2))+((3*h/8)*y(3))+((h/8)*y(4))-((3*(h^5)/80)*(y4(delta)));
end
y=1/pi*S;
end
can anyone fix it for me, to include the error... i got it right without it, but for some reason putting the error in, stuffs it up, and cant find where

it says something about the matrix must be square????

CB

4. here's the original qn

Let us consider an approximation to an integral. Let f(x) be some continuous function on
[a, b]. We wish to find an approximation for the integral
I = int from a to b of f(x)dx
in the following manner:
Subdivide the interval into N intervals of length h = (b−a)/N. Let xi = ih for i = 0, . . . ,N.
Let
Ij = int from 0 to h of f(xj+t)dt

Find a cubic polynomial Pj (x) that goes through (xj , f(xj)), (xj + h/3,f(xj + h/3)),(xj+2h/3, f(xj+2h/3) and (xj+1,f(xj+1))
We form an approximation for the integral by letting
I=sum(j=0 to N-1) of w0*f(xj)+w1*f(xj+h/3)+w2*f(xj+2h/3)+w3*f(xj+1)
Find these weights, wi.
In 2 peices of code, plot the first three Bessel functions, J0(x), J1(x) and J2(x), on the
interval [0, 20]. The first peice of code should be a MATLAB function BJ(x, n) outputing
the approximation for the integral representation of Jn, given by
Jn(x) =(1/pi)int from 0 to pi of cos(nt − x sin t)dt
using the above method for 100 subdivisions of [0, pi]. The second peice of code should call
the function an produce the required plots with 2000 subdivisions of [0, 20].

I dont need to find the error for qn, but i was just curious, on how you would do it? does it even matter?