Performing Newton's divided difference and trapezoidal rule in one program

How to run the code below without getting the error pictured below. The code is written like the Professor wants, but I cannot get past this error.

clear, close

%%%%NEWTON DIVIDED DIFFERENCE%%%

%all numbers in (cm)

l=50; %length

x=[2 4 6 8 10 12 14 16 8]; %x diameter in (cm)

fx=[5 7 8.4 9.1 9.1 8.3 8.1 7.8 6.6];%y diameter in (cm)

r=7; %radius of final bar

%INPUT Newtons divided difference

k=length(x);

B=zeros(k,k);

B(:,1)=fx;

for j=1:k-1

for i=1:k-j

B(i,j+1)=(B(i+1,j)-B(i,j))/(x(i+j)-x(i));

end

end

syms X

f=0;

m=1;

f=f+B(1,1)*m;

for i=1:k-1

m=m*(X-x(i));

f=f+B(1,i+1)*m;

end

p=3.3 ;

A=double(subs(f,X,p));

display(A)

syms U

f=0;

m=1;

f=f+B(1,1)*m;

for i=1:k-1

m=m*(U-x(i));

f=f+B(1,i+1)*m;

end

p=2 ;

A=double(subs(f,U,x));

display(A)

syms X;

f=B(1);

for i=1:k-1

f=f+B(i+1)*X^i; %f=required polynomial

end

disp(f);

syms X x

f=subs(f,X,x);

%Trapezodial Rule

%INPUT SECTION

f=@(x)f-(0.5*(pi*r^2)); %function

x=[2 4 6 8 10 12 14 16 8];

a=x(:,1); %lower limit

b=x(:,k); %upper limit

tol=1/100; %tolerance= percentage given/100

I0=0;

RE=100;

%CALCULATIONS

n=0;

while RE>tol

n=n+1;

h=(b-a)/n;

x_vec=a:h:b;

for i=1:length(x_vec)

fx(i)=f(x_vec(i));

end

I=(h/2)*(fx(1)+2*sum(fx(2:n))+fx(n+1));

RE=abs((I-I0)/I)*100;

fprintf('n1=%d \t I=%0.4f\t RE=%0.2f \n',[n;I;RE])

I0=I;

end

Answer 1

[Torsten]

I don't know if this is what you want. Your code is very badly structured and hard to follow.

Shouldn't

x=[2 4 6 8 10 12 14 16 18]; %x diameter in (cm)

instead of

x=[2 4 6 8 10 12 14 16 8]; %x diameter in (cm)

at the two places in your code ?

clear, close

%%%%NEWTON DIVIDED DIFFERENCE%%%

%all numbers in (cm)
l=50;   %length
x=[2 4 6 8 10 12 14 16 8]; %x diameter in (cm)
fx=[5 7 8.4 9.1 9.1 8.3 8.1 7.8 6.6];%y diameter in (cm)
r=7;    %radius of final bar

%INPUT  Newtons divided difference
k=length(x);
B=zeros(k,k);
B(:,1)=fx;

for j=1:k-1
    for i=1:k-j
        B(i,j+1)=(B(i+1,j)-B(i,j))/(x(i+j)-x(i));
    end
end




syms X

f=0;
m=1;
f=f+B(1,1)*m;

for i=1:k-1
    m=m*(X-x(i));
    f=f+B(1,i+1)*m;    
end

p=3.3  ;
A=double(subs(f,X,p));
display(A)

A = Inf

syms U

f=0;
m=1;

f=f+B(1,1)*m;



for i=1:k-1
    m=m*(U-x(i));
    f=f+B(1,i+1)*m;
    
end

p=2  ;
A=double(subs(f,U,x));
display(A)

A = 1x9

NaN NaN NaN NaN NaN NaN NaN NaN NaN

<mw-icon class=""></mw-icon>

<mw-icon class=""></mw-icon>

syms X;
f=B(1);
for i=1:k-1
    f=f+B(i+1)*X^i;      %f=required polynomial
end
disp(f);

syms X x
f=subs(f,X,x);



%Trapezodial Rule

%INPUT SECTION
f = matlabFunction(f);
f = @(x)f(x)-0.5*pi*r^2;
%f=@(x)f-(0.5*(pi*r^2));                    %function
x=[2 4 6 8 10 12 14 16 8];
a=x(:,1);                                  %lower limit
b=x(:,k);                                  %upper limit
tol=1/100;                                 %tolerance= percentage given/100

I0=0;
RE=100;

%CALCULATIONS
n=0;
while RE>tol
    n=n+1;
    h=(b-a)/n;
    x_vec=a:h:b;
    for i=1:length(x_vec)
        fxx(i)=f(x_vec(i));
    end
    I=(h/2)*(fxx(1)+2*sum(fxx(2:n))+fxx(n+1));
    RE=abs((I-I0)/I)*100;
    fprintf('n1=%d \t I=%0.4f\t RE=%0.2f \n',[n;I;RE])
    I0=I;
end

n1=1 	 I=388586617.3859	 RE=100.00 
n1=2 	 I=204334302.7859	 RE=90.17 
n1=3 	 I=158189801.3859	 RE=29.17 
n1=4 	 I=140867800.8710	 RE=12.30 
n1=5 	 I=132632681.7456	 RE=6.21 
n1=6 	 I=128100325.1859	 RE=3.54 
n1=7 	 I=125347179.1992	 RE=2.20 
n1=8 	 I=123552066.5293	 RE=1.45 
n1=9 	 I=122317591.8700	 RE=1.01 
n1=10 	 I=121432702.3593	 RE=0.73 
n1=11 	 I=120776974.6596	 RE=0.54 
n1=12 	 I=120277664.3085	 RE=0.42 
n1=13 	 I=119888738.9331	 RE=0.32 
n1=14 	 I=119579923.8198	 RE=0.26 
n1=15 	 I=119330648.9816	 RE=0.21 
n1=16 	 I=119126543.6654	 RE=0.17 
n1=17 	 I=118957323.3025	 RE=0.14 
n1=18 	 I=118815471.0917	 RE=0.12 
n1=19 	 I=118695390.2871	 RE=0.10 
n1=20 	 I=118592844.5417	 RE=0.09 
n1=21 	 I=118504579.6678	 RE=0.07 
n1=22 	 I=118428062.6072	 RE=0.06 
n1=23 	 I=118361297.8758	 RE=0.06 
n1=24 	 I=118302696.2804	 RE=0.05 
n1=25 	 I=118250979.5659	 RE=0.04 
n1=26 	 I=118205110.1803	 RE=0.04 
n1=27 	 I=118164238.8727	 RE=0.03 
n1=28 	 I=118127665.1305	 RE=0.03 
n1=29 	 I=118094806.9846	 RE=0.03 
n1=30 	 I=118065177.7305	 RE=0.03 
n1=31 	 I=118038367.8118	 RE=0.02 
n1=32 	 I=118014030.5997	 RE=0.02 
n1=33 	 I=117991871.1385	 RE=0.02 
n1=34 	 I=117971637.1709	 RE=0.02 
n1=35 	 I=117953111.9307	 RE=0.02 
n1=36 	 I=117936108.3159	 RE=0.01 
n1=37 	 I=117920464.1472	 RE=0.01 
n1=38 	 I=117906038.2886	 RE=0.01 
n1=39 	 I=117892707.4542	 RE=0.01 
n1=40 	 I=117880363.5659	 RE=0.01 
n1=41 	 I=117868911.5571	 RE=0.01