Plotting from For Loop ODEFUN and BVP4C Initial Guess

3 次查看（过去 30 天）
Tony Stianchie
Tony Stianchie 2023-4-29
评论： Torsten 2023-4-29
%% 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.
  4 个评论
显示 2更早的评论
Star Strider
Star Strider 2023-4-29
Ran in:
‘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
Tony Stianchie 2023-4-29
%% 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

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Torsten
Torsten 2023-4-29
Ran in:
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
  4 个评论
显示 2更早的评论
Torsten
Torsten 2023-4-29
Ran in:
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

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Boundary Value Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by