# what is wrong with this code?

2 views (last 30 days)
Vetrichelvan Pugazendi on 18 Apr 2018
Commented: Walter Roberson on 25 Apr 2018
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

Veera Kanmani on 18 Apr 2018
what is the error here? Can you please specify the error?
Vetrichelvan Pugazendi on 18 Apr 2018
The error is undefined variable ‘int’ in the ending

Walter Roberson on 18 Apr 2018

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 on 25 Apr 2018
Edited: Walter Roberson on 25 Apr 2018

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 on 25 Apr 2018
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?