Error using double (subs())

9 views (last 30 days)
rahul rai
rahul rai on 3 Aug 2022
Answered: Torsten on 3 Aug 2022
Hi,
I want to run newton's method and I am not able to convert the syms variable into double even though I am using subs.
Here is my code.
syms initial_R
u = -100;
v = @(initial_R) initial_R.*exp ((2*(1-alpha))/(2+(1-alpha)*u.*initial_R)./(2-alpha*u.*initial_R).^((alpha)^2) * (2+(1-alpha)*u.*initial_R).^(1-alpha^2));
vpr = diff (v,initial_R);
funcr = @(initial_R) vpr;
kmax = 10;
initial_R(1) = 0;
initial_R(2) = 0.0070;
y(2) = v(initial_R(2));
ypr(2) = funcr (initial_R(2));
tol = 0.000001
for k = 2:0.01:kmax
initial_R(k+1) = initial_R(k) - y(k)/ypr(k);
y(k+1) = feval (v,initial_R(k+1));
c = initial_R(k+1) - initial_R(k);
d = subs(c, initial_R(k+1), initial_R(k));
if (abs(double(d))) < tol
disp('The Numerical has ended')
end
ypr(k+1) = feval(funcr, initial_R(k+1));
end
I would like to know how else can I check the tolerance level because everytime I run this code with no matter how many variations, the if line throws an error.
Error using symengine
Unable to convert expression containing symbolic variables into double array. Apply 'subs' function first to substitute values for
variables.
Error in sym/double (line 868)
Xstr = mupadmex('symobj::double', S.s, 0);
Related documentation

Answers (1)

Torsten
Torsten on 3 Aug 2022
syms initial_R
u = -100;
alpha = 2;
v = initial_R*exp((2*(1-alpha))/(2+(1-alpha)*u*initial_R)/(2-alpha*u*initial_R)^(alpha)^2*(2+(1-alpha)*u*initial_R)^(1-alpha^2))
v = 
vpr = diff(v,initial_R);
y = matlabFunction(v)
y = function_handle with value:
@(initial_R)initial_R.*exp(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^4.*-2.0)
ypr = matlabFunction(vpr)
ypr = function_handle with value:
@(initial_R)exp(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^4.*-2.0)+initial_R.*exp(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^4.*-2.0).*(1.0./(initial_R.*1.0e+2+2.0).^4.*1.0./(initial_R.*2.0e+2+2.0).^5.*1.6e+3+1.0./(initial_R.*1.0e+2+2.0).^5.*1.0./(initial_R.*2.0e+2+2.0).^4.*8.0e+2)
kmax = 50;
initial_R_old = 1.0;
tol = 0.000001;
k = 0;
dx = 1.0;
while dx > tol && k < kmax
k = k+1;
initial_R_new = initial_R_old - y(initial_R_old)/ypr(initial_R_old);
dx = abs(initial_R_old - initial_R_new);
initial_R_old = initial_R_new;
end
initial_R_new
initial_R_new = 0
y(initial_R_new)
ans = 0

Products


Release

R2022a

Community Treasure Hunt

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

Start Hunting!