what is wrong with this code?

I have been working on this code and i encountered this unsolvable error. Please help !

code :

clc

clear all 

close all 

%%Parameters of Blasius Equation

U_inf = 1;

L = 10;

mu = 1.789E-5;

rho = 1.225;

nu = mu/rho;

A = sqrt(nu/U_inf);

h = 0.01;

%%Numerical Solution of Blasius Equation Using Runge-Kutta

f1 = @(x, y1, y2, y3) y2;

f2 = @(x, y1, y2, y3) y3;

f3 = @(x, y1, y2, y3) -y1*y3;

eta = 0:h:10;

x = 0:h:10;

y1(1) = 0;

y2(1) = 0;

y3(1) = 0.4696;

for i = 1:(length(eta)-1)

a = h.*[f1(eta(i), y1(i), y2(i), y3(i)), f2(eta(i), y1(i), y2(i), y3(i)), f3(eta(i), y1(i), y2(i), y3(i))];

b = h.*[f1(eta(i), y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f2(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f3(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2)];

c = h.*[f1(eta(i), y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f2(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f3(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2)];

d = h.*[f1(eta(i), y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f2(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f3(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3))];

y3(i+1) = y3(i)+ 1/6*(a(3)+2*b(3)+2*c(3)+d(3));

y2(i+1) = y2(i)+ 1/6*(a(2)+2*b(2)+2*c(2)+d(2));

y1(i+1) = y1(i)+ 1/6*(a(1)+2*b(1)+2*c(1)+d(1));

end

%%Plotting and Visualization

figure(1)

plot(eta,y1,eta, y2, eta, y3, 'LineWidth', 2)

xlim([0 10])

title('Solution of Blasius eqution', 'FontSize', 14);

xlabel('f, f'' and f''''', 'FontSize', 20);

ylabel('\eta', 'FontSize', 20);

grid on

Legend1 = {'f(\eta)', 'f''(\eta)', 'f''''(\eta)'};

legend(Legend1, 'FontSize', 14);

syms q a b pr a1 b1;
pr=1;
a=((diff([y2],2)).^(pr));
b=((diff([y2],2)).^(pr));
a1=int(a,'eta',0,10);
b1=int(b,'eta',0,inf);
q=a1/b1

Answer 1

[Veera Kanmani]

what is the error here? Can you please specify the error?

[Vetrichelvan Pugazendi]

The error is undefined variable ‘int’ in the ending

Answer 2

[Walter Roberson]

Your y2 is numeric, so diff(y2, 2) is requesting the second numeric differences of the numeric array. The resulting numeric array is taken to the power of the numeric variable pr (value 1) and the result is stored in a. Then b is calculated exactly the same way, giving exactly the same result.

So now a is a pure numeric array. And you seem to intend to integrate it with respect to the symbolic variable eta. But there is no routine named int() that accepts numeric first argument.

What you could do is to

int(sym(a), 'eta', 0, 10)

There would be little point in doing that, but it could be done. The integral of a numeric constant with respect to definite bounds is just the constant times the difference between the bounds, so just multiplying a by 10 would achieve the same effect.

[Vetrichelvan Pugazendi]

I changed the coding to this and tried it

Code :

clc

clear all

close all

%% Parameters of Blasius Equation

U_inf = 1;

L = 10;

mu = 1.789E-5;

rho = 1.225;

nu = mu/rho;

A = sqrt(nu/U_inf);

h = 0.01;

%% Numerical Solution of Blasius Equation Using Runge-Kutta

f1 = @(x, y1, y2, y3) y2;

f2 = @(x, y1, y2, y3) y3;

f3 = @(x, y1, y2, y3) -y1*y3;

eta = 0:h:10;

x = 0:h:10;

y1(1) = 0;

y2(1) = 0;

y3(1) = 0.4696;

for i = 1:(length(eta)-1)

a = h.*[f1(eta(i), y1(i), y2(i), y3(i)), f2(eta(i), y1(i), y2(i), y3(i)), f3(eta(i), y1(i), y2(i), y3(i))];

b = h.*[f1(eta(i), y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f2(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f3(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2)];

c = h.*[f1(eta(i), y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f2(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f3(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2)];

d = h.*[f1(eta(i), y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f2(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f3(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3))];

y3(i+1) = y3(i)+ 1/6*(a(3)+2*b(3)+2*c(3)+d(3));

y2(i+1) = y2(i)+ 1/6*(a(2)+2*b(2)+2*c(2)+d(2));

y1(i+1) = y1(i)+ 1/6*(a(1)+2*b(1)+2*c(1)+d(1));

end

%% Plotting and Visualization

figure(1)

plot(eta,y1,eta, y2, eta, y3, 'LineWidth', 2)

xlim([0 10])

title('Solution of Blasius eqution', 'FontSize', 14);

xlabel('f, f'' and f''''', 'FontSize', 20);

ylabel('\eta', 'FontSize', 20);

grid on

Legend1 = {'f(\eta)', 'f''(\eta)', 'f''''(\eta)'};

legend(Legend1, 'FontSize', 14);

syms q pr a1 b1; 
pr=1; 
%a=((diff([y2],2)).^(pr)); 
%b=((diff([y2],2)).^(pr)); 
a1= int(sym(((diff([y2],2)).^(pr))), 'eta', 0, 10);
b1= int(sym(((diff([y2],2)).^(pr))), 'eta' , 0, inf);
q=a1/b1

Still there’s an error

Please try copy pasting the code and running it

[Walter Roberson]

The integral of a numeric constant with respect to definite bounds is just the constant times the difference between the bounds. Your sym(((diff([y2],2)).^(pr)) is a vector of numeric constants. Your bounds are 0 and inf. Numeric constant times difference in bounds is numeric constant times (inf-0) which is numeric constant times inf which is inf for all positive constants and -inf for all negative constants, and 0 if the constant is 0.

It happens that your y2 includes negative values and some zeros. Therefore your b1 consists entirely of -inf and 0.

Your q=a1/b1 is a matrix right division operation, roughly equivalent to

q = a1 * pinv(b1)

However the vector of -inf and 0 mixed is calculated as being inconsistent with a1 so you get a warning message displayed about inconsistency. Then after that the symbolic toolbox runs into an internal error in how it deals with the inconsistent system.

It is pretty unlikely that you should be using int() with those numeric constants. I suspect you should be using trapz() to do a numeric integration.

[Vetrichelvan Pugazendi]

I still have an error saying - Error using sym/pinv, argument must not contain ‘Nan’ or ‘inf’

[Walter Roberson]

What result would you expect for int(sym(2), 'eta', 0, inf) ?

[Vetrichelvan Pugazendi]

Fix the error please

[Walter Roberson]

In order fix the error you need to tell me what your intended result is for the line where you wrote int() of a numeric value from 0 to infinity. Suppose your b value came out as 2, then what would you have liked the integral to be? The integral of any nonzero constant over an infinite range is always infinite.

It might help to post an image of the formula that you are trying to implement.

[Vetrichelvan Pugazendi]

I would like you to solve for energy equation and i have solved for momentum equation, using momentum equation solve for energy equation.

This is the code for momentum equations.

clc

clear all 

close all 

U_inf = 1;

L = 10;

mu = 1.789E-5;

rho = 1.225;

nu = mu/rho;

A = sqrt(nu/U_inf);

h = 0.01;

%%  Runge-Kutta Method

f1 = @(x, y1, y2, y3) y2;

f2 = @(x, y1, y2, y3) y3;

f3 = @(x, y1, y2, y3) -y1*y3;

eta = 0:h:10;

x = 0:h:10;

y1(1) = 0;

y2(1) = 0;

y3(1) = 0.4696;

for i = 1:(length(eta)-1)

a = h.*[f1(eta(i), y1(i), y2(i), y3(i)), f2(eta(i), y1(i), y2(i), y3(i)), f3(eta(i), y1(i), y2(i), y3(i))];

b = h.*[f1(eta(i), y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f2(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2), f3(eta(i)+h/2, y1(i)+a(1)/2, y2(i)+a(2)/2, y3(i)+a(3)/2)];

c = h.*[f1(eta(i), y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f2(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2), f3(eta(i)+h/2, y1(i)+b(1)/2, y2(i)+b(2)/2, y3(i)+b(3)/2)];

d = h.*[f1(eta(i), y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f2(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3)), f3(eta(i)+h, y1(i)+c(1), y2(i)+c(2), y3(i)+c(3))];

y3(i+1) = y3(i)+ 1/6*(a(3)+2*b(3)+2*c(3)+d(3));

y2(i+1) = y2(i)+ 1/6*(a(2)+2*b(2)+2*c(2)+d(2));

y1(i+1) = y1(i)+ 1/6*(a(1)+2*b(1)+2*c(1)+d(1));

end

%% Graph Result

figure(1)

plot(eta,y1,eta,y2,eta,y3, 'LineWidth', 2)

xlim([0 10])

title('Solution of Blasius eqution', 'FontSize', 14);

ylabel('f, f'' and f''''', 'FontSize', 20);

xlabel('\eta', 'FontSize', 20);

grid on

Legend1 = {'f(\eta)', 'f''(\eta)', 'f''''(\eta)'};

legend(Legend1, 'FontSize', 14);

[Walter Roberson]

Those look like differential equations to me. Numeric solutions for differential equations are typically most easily solved using one of the ode*() routines such as ode45().

I do not seem to recognize the meaning of "Pr." in the energy equation?