newton raphson methon and symbolic functions
1 Ansicht (letzte 30 Tage)
Ältere Kommentare anzeigen
Marko
am 28 Dez. 2023
Beantwortet: Sulaymon Eshkabilov
am 28 Dez. 2023
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 Kommentare
Akzeptierte Antwort
Torsten
am 28 Dez. 2023
Bearbeitet: Torsten
am 28 Dez. 2023
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 Kommentare
Weitere Antworten (2)
Hassaan
am 28 Dez. 2023
% 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 Kommentare
Sulaymon Eshkabilov
am 28 Dez. 2023
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 Kommentare
Siehe auch
Kategorien
Mehr zu Symbolic Math Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!