Simpsons 3.8 rule, error help matlab

May 2009
45
1
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????
 
Last edited by a moderator:
May 2009
45
1
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?
 

CaptainBlack

MHF Hall of Fame
Nov 2005
14,972
5,271
someplace
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????
Please post the original question.

CB
 
  • Like
Reactions: fredrick08
May 2009
45
1
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?
 
Similar Math Discussions Math Forum Date
Calculus
Calculus
Calculus
Calculus