Plotting from For Loop ODEFUN and BVP4C Initial Guess

%% Thetadot(0) vs Pr 

m = 0;
Prinf = 1000;

for Pr = linspace(0,Prinf,1)

solinit = bvpinit(Pr,[0 0 0 0 0.05]);
sol = bvp4c(@odefun, @odefun_bc,solinit);

xint = linspace(0,Prinf,1);
Sxint = deval(sol,xint);

figure(19)
hold on
title('HeatFlux(0) vs Pr')
xlabel('Pr')
ylabel('Heat Flux')
plot(xint,Sxint(5,1)); % plots qdot(0)
end

I'd like to plot y(5,1) across varying Pr; however it seems my initial guess is not sufficient.

[Torsten]

Is this really what you want ?

Prinf = 1000;
linspace(0,Prinf,1)

ans = 1000

[Tony Stianchie]

No, that's a typo. I'd like to generate a vector (0 to 1000) with n =1 incremenets

[Star Strider]

‘I'd like to generate a vector (0 to 1000) with n =1 incremenets’

Prinf = 1000;
xint = linspace(0,Prinf,Prinf+1)

xint = 1×1001

0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29

.

[Tony Stianchie]

%% Thetadot(0) vs Pr 

m = 1;
Prinf = 1000;
x = zeros(1,Prinf);
   
for Pr = linspace(0,Prinf,Prinf+1)

solinit = bvpinit(linspace(0,etainf,100),[0 0 0 1 0]);
sol = bvp4c(@odefun, @odefun_bc,solinit);

xint = linspace(0,etainf,100);
Sxint = deval(sol,xint);

x(1,Pr+1) = Sxint(5,1);

end

figure(3)
plot(linspace(0,Prinf,Prinf+1), x)

Thanks - I wrote the for loop as such. I'm using the same odefun as attached previously

Answer 1

[Torsten]

global m Pr 

etainf = 20; % Find Convergence for both Temp and Velocity

%% Thetadot(0) vs Pr 

m = 0;
PR = 0:0.5:50;
for i = 1:numel(PR)
  Pr = PR(i);
  solinit = bvpinit(linspace(0,etainf,100),[0 0 0 0 0.05]);
  sol = bvp4c(@odefun, @odefun_bc,solinit);
  qdot0(i) = sol.y(5,1);
end
plot(PR,qdot0)
title('HeatFlux(0) vs Pr')
xlabel('Pr')
ylabel('Heat Flux')

function yprime = odefun(eta,y)

yprime = zeros(5,1);

global m Pr

% Blasius Eqn
yprime(1) = y(2);
yprime(2) = y(3);
yprime(3) = (-1/2)*(m+1)*y(1)*y(3)+m*y(2)^2 - m;

% Energy Eqn
yprime(4) = y(5);
yprime(5) = (-1/2)*Pr*y(1)*(m+1)*y(5);


yprime = yprime';
end

function res = odefun_bc(ya, yb) 
res = [ya(1); ya(2); yb(2)-1; ya(4); yb(4)-1];
end

[Tony Stianchie]

If I wanted to fit a curve to the plotted graph, how would you recommend I proceed?

[Torsten]

And what is (are) the fitting parameter(s) ?

[Tony Stianchie]

I expect it to plot to a polynomial

[Torsten]

Looks like a + b*sqrt(Pr) in my opinion.

global m Pr 
etainf = 20; % Find Convergence for both Temp and Velocity
%% Thetadot(0) vs Pr 
m = 0;
PR = 0:0.5:50;
for i = 1:numel(PR)
    Pr = PR(i);
    solinit = bvpinit(linspace(0,etainf,100),[0 0 0 0 0.05]);
    sol = bvp4c(@odefun, @odefun_bc,solinit);
    qdot0(i) = sol.y(5,1);
end
hold on
plot(PR,qdot0)
A = [ones(numel(PR),1) sqrt(PR.')];
b = qdot0.';
x = A\b;
plot(PR,x(1)+x(2)*sqrt(PR))
hold off

function yprime = odefun(eta,y)
yprime = zeros(5,1);
global m Pr
% Blasius Eqn
yprime(1) = y(2);
yprime(2) = y(3);
yprime(3) = (-1/2)*(m+1)*y(1)*y(3)+m*y(2)^2 - m;
% Energy Eqn
yprime(4) = y(5);
yprime(5) = (-1/2)*Pr*y(1)*(m+1)*y(5);
yprime = yprime';
end
function res = odefun_bc(ya, yb) 
res = [ya(1); ya(2); yb(2)-1; ya(4); yb(4)-1];
end