Plotting from For Loop ODEFUN and BVP4C Initial Guess

2023 4 29
1 个回答
回答已采纳
更新时间：2023 4 29
3 次查看（30 天）

0 个投票

%% 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更早的评论 隐藏 2更早的评论

Torsten
Torsten 2023-4-29
Ran in:
Ran in:
Is this really what you want ?
Prinf = 1000;
linspace(0,Prinf,1)
ans = 1000
Tony Stianchie
Tony Stianchie 2023-4-29
No, that's a typo. I'd like to generate a vector (0 to 1000) with n =1 incremenets
Star Strider
Star Strider 2023-4-29
Ran in:
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:

1 个投票

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更早的评论 隐藏 2更早的评论

Tony Stianchie
Tony Stianchie 2023-4-29
If I wanted to fit a curve to the plotted graph, how would you recommend I proceed?
Torsten
Torsten 2023-4-29
编辑：Torsten 2023-4-29
And what is (are) the fitting parameter(s) ?
Tony Stianchie
Tony Stianchie 2023-4-29
I expect it to plot to a polynomial
Torsten
Torsten 2023-4-29
Ran in:
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 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Simscape Multibody 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by