Azzera filtri
Azzera filtri

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

2 visualizzazioni (ultimi 30 giorni)
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

Risposta accettata

Torsten
Torsten il 15 Lug 2024
Modificato: Torsten il 15 Lug 2024
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

Più risposte (0)

Categorie

Scopri di più su Symbolic Math Toolbox in Help Center e File Exchange

Prodotti


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by