newton raphson methon and symbolic functions
์กฐํ ์: 1 (์ต๊ทผ 30์ผ)
์ด์ ๋๊ธ ํ์
Marko
2023๋
12์ 28์ผ
๋ต๋ณ: Sulaymon Eshkabilov
2023๋
12์ 28์ผ
Hello, I have to write a code that calculates how much the spherical tank must be filled to contain 30m3 if R=3m using Newton Raphson method. The volume is calculated according to the expression ๐=๐*โ^2* [3*๐
โโ]/ 3. I have to print the aproximations and relative error. I have tried writing the code multipule times, but Ialways have errors with derivative of a function and/or solving the equation for the exact answer (functions diff and solve). This is my code so far:
R=3;
h=0:0.1:6;
V=pi.*h.^2.*(3*R-h)/3;
plot(h,V,'k');
V_desired=30;
x=3;
max_iterations=100;
tolerance=1e-6;
answ=solve('pi*h^2*(3*3-h/3==V_desired', h);
A(1,1)=('iteration');
A(1,2)=('aproximation');
A(1,3)=('relative error');
for iteration=1:max_iteration
H=volume(x, R);
HD=dvolume(x,R);
x=x-H/HD;
error=(answ-x)/answ*100;
A(iteration+1, 1)=iteration;
A(iteration+1, 2)=x;
A(iteration+1, 3)=error;
if error<tolerance
break
end
end
disp(A);
This is a file for function:
function vol=volume(h,R)
vol=pi*h^2*(3*R-h)/3;
end
And this one is for derivative of that function:
functiondvol=dvolumen(h,R)
dvol=diff(volume(h,R), h)
end
Any answer would help, thank you!
๋๊ธ ์: 0
์ฑํ๋ ๋ต๋ณ
Torsten
2023๋
12์ 28์ผ
ํธ์ง: Torsten
2023๋
12์ 28์ผ
Better differentiate your function with respect to h in advance and pass the result to "dvolume" so that you only need to substitute the actual h value.
Further I suggest using a while-loop instead of a for-loop and to add the computation of the error in volume(h), not only in h.
R=3;
h=0:0.1:6;
V=pi.*h.^2.*(3*R-h)/3;
plot(h,V,'k');
V_desired=30;
xold=3;
max_iterations=100;
tolerance=1e-6;
%answ=solve('pi*h^2*(3*R-h)/3==V_desired', h);
%A(1,1)=('iteration');
%A(1,2)=('aproximation');
%A(1,3)=('relative error');
A(1,1) = 0;
A(1,2) = xold;
A(1,3) = NaN;
for iteration=1:max_iterations
H=volume(xold, R, V_desired);
HD=dvolume(xold,R, V_desired);
x=xold-H/HD;
error=abs((xold-x)/xold)*100;
A(iteration+1, 1)=iteration;
A(iteration+1, 2)=x;
A(iteration+1, 3)=error;
xold = x;
if error<tolerance
break
end
end
disp(A);
This is a file for function:
function vol=volume(h,R, V_desired)
vol=pi*h^2*(3*R-h)/3-V_desired;
end
And this one is for derivative of that function:
function dvol=dvolume(h,R, V_desired)
%dvol=diff(volume(h,R), h)
syms hsym
dvol = diff(pi*hsym^2*(3*R-hsym)/3-V_desired,hsym);
dvol = double(subs(dvol,hsym,h));
end
๋๊ธ ์: 0
์ถ๊ฐ ๋ต๋ณ (2๊ฐ)
Hassaan
2023๋
12์ 28์ผ
% Define the radius of the tank and the desired volume
R = 3;
v_desired = 30;
% Initial guess for h
h = 1;
% Define tolerance and maximum number of iterations
tolerance = 1e-6;
max_iterations = 100;
% Newton-Raphson iteration
for iteration = 1:max_iterations
current_volume = volume(h, R);
current_derivative = dvolume(h, R);
h_new = h - (current_volume - v_desired) / current_derivative;
% Calculate relative error
error = abs((h_new - h) / h_new);
% Display current iteration, h value, and error
fprintf('Iteration %d: h = %.6f, Error = %.6f\n', iteration, h_new, error);
% Check if the error is within the tolerance
if error < tolerance
break;
end
% Update h for the next iteration
h = h_new;
end
% Display the final height
fprintf('\nFinal height h = %.6f after %d iterations\n', h, iteration);
% Volume function
function vol = volume(h, R)
vol = pi * h^2 * (3 * R - h) / 3;
end
% Corrected derivative function
function dvol = dvolume(h, R)
syms hs % Define a symbolic variable for h
vol = pi * hs^2 * (3 * R - hs) / 3; % Define the volume as a symbolic expression
dvolSym = diff(vol, hs); % Take the symbolic derivative with respect to hs
dvol = subs(dvolSym, hs, h); % Substitute the symbolic variable with the numeric value of h
dvol = double(dvol); % Convert symbolic result to numeric
end
------------------------------------------------------------------------------------------------------------------------------------------------
If you find the solution helpful and it resolves your issue, it would be greatly appreciated if you could accept the answer. Also, leaving an upvote and a comment are also wonderful ways to provide feedback.
๋๊ธ ์: 0
Sulaymon Eshkabilov
2023๋
12์ 28์ผ
Here is another alt. solution:
R=3;
h=0:0.1:6;
V=pi.*h.^2.*(3*R-h)/3;
plot(h,V,'k');
IDX = find(ceil(V)==30);
hold on
plot(h(IDX), V(IDX), 'rd', 'MarkerSize', 13, 'MarkerFaceColor','m')
legend('V', 'Solution point of h')
V_desired=30;
x=0.5; % Initial guess value
max_iterations=100;
tolerance=1e-6;
syms hh
EQN = pi*hh^2*(3*R-hh)/3==V_desired;
answ=solve(EQN, hh);
Answer = answ(2);
Answer = double(Answer);
disp(['Iters ', ' Appr ', ' Error '])
for iteration=1:max_iterations
Error=(Answer-x)/Answer;
A(iteration, :)=[iteration, x, abs(Error)];
fprintf('%d %15.6f %15.7f \n', A(iteration,:));
[Vol, dVol]=VOLUME(x);
x=x-Vol/dVol;
if abs(Error)<tolerance
disp('Iteration is halted becasue error tolerance is reached')
break
end
end
function [Vol, dVol]=VOLUME(h)
R = 3;
V_desired=30;
Vol=pi*h^2*(3*R-h)/3-V_desired;
syms hh
dV=diff(pi*hh^2*(3*R-hh)/3-30-V_desired, hh);
dVol=double(subs(dV, hh,h));
end
๋๊ธ ์: 0
์ฐธ๊ณ ํญ๋ชฉ
์นดํ ๊ณ ๋ฆฌ
Help Center ๋ฐ File Exchange์์ Symbolic Math Toolbox์ ๋ํด ์์ธํ ์์๋ณด๊ธฐ
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!