Method of Characteristics for Nozzle: Error with syms function

I am trying to create a Method of characteristics Code from JD Andersons book "Modern Compressible Flow with Historical Perspective" using example 11.1. Error arises when I try to obtain the coordinates of the points 2-7; where 7 is the n_div which I have taken.

Code:

clear;

clc;

close all;

%% Initial parameters %%

radian_to_deg = 180/pi;

degree_to_radian = pi/180;

prompt = "Enter a Mach Number ranging from 1.5 to 4: ";

Me = input(prompt);

prompt_1 = "Enter the number of divisions: ";

n_div = input(prompt_1);

%prompt_2 = "Enter the throat radius: ";

%TR = input(prompt_2);

g = 1.4;

mach_angle= [];

theta = [];

theta_a = [];

y=[];

x=[];

v=[];

v_a=[];

dydx_minus = [];

dydx_plus = [];

new_Mach=[];

%% get mach angle

ma_f = @(x) (asind(1/x)) ;

%% get theta_max

A = sqrt((g+1)/(g-1));

B = (g-1)/(g+1);

pm_f = @(x) (A*atand(sqrt(B*(x^2-1))) - atand(sqrt(x^2-1)));

theta_max = pm_f(Me)/2;

%% extraction of fractional part of theta_max

theta_i = sign(theta_max)*(abs(theta_max) - floor(abs(theta_max)));

%% Storing different values of theta_a

for i = 1:n_div

theta_a(i,1) = theta_i ;

theta_i = theta_i + 3.0;

end

v_a = theta_a;

%% Getting properties of a-1 characteristic line

ya=1;

xa=0;

%theta_a(1,1) = theta_i(1,1);

v_a(1,1) = theta_a(1,1);

%point 1

y(1,1) = 0;

Km_1 = theta_a(1,1) + v_a(1,1);

theta(1,1)= 0.5*(Km_1);

v(1,1) = 0.5*(Km_1);

%% get mach number for v(1,1)

syms x

f = @(x) (A*atand(sqrt(B*(x^2-1))) - atand(sqrt(x^2-1))-v(1,1));

f1 = matlabFunction(diff(f(x)));

%NR Method using a while loop

change = 10;

Mguess = 1.1;

while (change>1e-6)

Mnew = Mguess - f(Mguess)/f1(Mguess);

change = abs(Mguess - Mnew);

Mguess = Mnew;

end

new_Mach(1,1) = Mnew;

mach_angle(1,1) = ma_f(Mnew);

x(1,1) = -ya/tand(theta(1,1)-mach_angle(1,1));

%% rest of points:

for i = 2:n_div

Km = (theta_a(i,1)) + (v_a(i,1));

Kp = theta(i-1,1) - v(i-1,1);

theta(i,1) = 0.5*(Km+Kp);

v(i,1) = 0.5*(Km-Kp);

end

%% Calculating the Unknown Mach Number and corresponding Mach Angle using NR Method

syms x

for i=2:n_div

f = @(x) (A*atand(sqrt(B*(x^2-1))) - atand(sqrt(x^2-1))-v(i,1));

f1 = matlabFunction(diff(f(x)));

Mguess = 1.1;

change = 10;

while (change>1e-6)

Mnew = Mguess - f(Mguess)/f1(Mguess);

change = abs(Mguess - Mnew);

Mguess = Mnew;

end

new_Mach(i,1) = Mnew;

mach_angle(i,1) = ma_f(Mnew);

end

%% slopes and locations of points from 2-7

for i=2:n_div

dydx_minus(i,1) = tand(0.5*(theta(i)+theta(i)) - 0.5*(mach_angle(i)+mach_angle(i)));

dydx_plus(i,1) = tand(0.5*(theta(i) + theta(i-1))) + (0.5*(mach_angle(i) + mach_angle(i-1)));

x(i,1) = (y(i-1,1) - ya - (dydx_plus(i,1)*x(i-1)))/(dydx_minus(i,1) - dydx_plus(i,1));

y(i,1) = (dydx_plus(i,1)*(x(i,1)-x(i-1))) + y(i-1);

end

%slope and location of wall point 8

slopew1_minus = tand(theta_max);

slopew1_plus = dydx_plus(end,1);

xw = (y(end,1)-ya-(slopew1_plus*x(end,1)))/(slopew1_minus - slopew1_plus);

yw = (slopew1_plus*(xw(1,1)-x(end,1)))+y(end,1);

Error:

Error in MOC_prac (line 111)

y(i,1) = (dydx_plus(i,1)*(a_double(i,1)-a_double(i-1))) + y(i-1);

Unable to perform assignment because value of type 'sym' is not convertible to 'double'.

Error in MOC_prac (line 111)

y(i,1) = (dydx_plus(i,1)*(x(i,1)-x(i-1))) + y(i-1);

Caused by:

Error using mupadengine/feval2char

Unable to convert expression containing symbolic variables into double array. Apply 'subs' function

first to substitute values for variables.

Is there any way I can resolve this error? Thank you!

Answer 1

[Steven Lord]

Looking at a few sections in your code, first you define x to be a symbolic variable and you use it as such.

syms x
f = @(x) (A*atand(sqrt(B*(x^2-1))) - atand(sqrt(x^2-1))-v(1,1));
f1 = matlabFunction(diff(f(x)));

Then about 10 lines later you assign into it. Since that assignment is indexed assignment the value on the right side is converted into the sym class and stored in x, so x remains a sym array.

x(1,1) = -ya/tand(theta(1,1)-mach_angle(1,1));

Seven lines later you again state that x is a symbolic variable and use it as such.

syms x
for i=2:n_div
f = @(x) (A*atand(sqrt(B*(x^2-1))) - atand(sqrt(x^2-1))-v(i,1));
f1 = matlabFunction(diff(f(x)));

Then you store a value to that symbolic variable using indexed assignment again, which means x stays a sym object. One point to note is that nowhere here do you assign a value to x(1) so that will remain a plain old 'x'.

for i=2:n_div
    dydx_minus(i,1) = tand(0.5*(theta(i)+theta(i)) - 0.5*(mach_angle(i)+mach_angle(i)));
    dydx_plus(i,1) = tand(0.5*(theta(i) + theta(i-1))) + (0.5*(mach_angle(i) + mach_angle(i-1)));
    x(i,1) = (y(i-1,1) - ya - (dydx_plus(i,1)*x(i-1)))/(dydx_minus(i,1) - dydx_plus(i,1));

Now you use that symbolic variable in the computations on the right-hand side. This means that right hand side is symbolic. MATLAB will need to try to convert that into double (because you're performing indexed assignment into y on the left hand side) but it can't.

    y(i,1) = (dydx_plus(i,1)*(x(i,1)-x(i-1))) + y(i-1);

Because x(1) is just 'x', the right-hand side is likely to be something like <number>*x + <number>. Without a value for x, this cannot be converted to a number.

If I asked you to tell me what number 2*x + 1 is, your first question is likely to be "What is x?" MATLAB can't ask you "What is x?" here so it just says essentially "I can't tell you what number 2*x+1 is."

Use subs to substitute in a value for the symbolic variable x in the x vector [you may want to use a different variable name for the variable you use as x(i, 1) than the one you use as x(i-1)] so MATLAB can convert that expression into double.