# Method of Characteristics for Nozzle: Error with syms function

23 ビュー (過去 30 日間)
Bamelari Jovani 2024 年 7 月 18 日 13:57

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 %%
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:
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!

サインインしてコメントする。

### 採用された回答

Steven Lord 2024 年 7 月 18 日 15:10
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.

サインインしてコメントする。

### カテゴリ

Help Center および File ExchangeConversion Between Symbolic and Numeric についてさらに検索

R2024a

### Community Treasure Hunt

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

Start Hunting!

Translated by