- Initialize Data Storage: Store the results for each value of Prf in a matrix.
- Loop Over Prf: Solve the boundary value problem for each value of Prf.
- Create the Surface Plot: Use surf to plot the data.
Please help me to run surf plot with bvp4c.
34 次查看(过去 30 天)
显示 更早的评论
Please help me to run surf plot with bvp4c.The surfce digram consists of (constant Prf [2 :1:6] represents y-axis & vector>> sol.x [0 4] represents x-axis (m = linspace(0,4);) & second solution only sol.y(6,:) represents z-axis).The following is the code for 2D (sol.x [0 4] and only sol.y(6,:)). How to give command for making surf plot in bvp4c.
proj()
function sol= proj
clc;clf;clear;
%Relation of base fluid
rhof=1;kf=0.613*10^5;cpf=4179*10^4;muf=10^-3*10;sigf=0.05*10^-8;alfaf=kf/(rhof*cpf);
%FE3O4
ph1=0.01;rho1=5180*10^-3;cp1=670*10^4;k1=9.7*10^5;sig1=0.74*10^-2;
%copper
ph2=0.01;rho2=8933*10^-3;cp2=385*10^4;k2=401*10^5;sig2=5.96*10^-1;
%Relation of hyprid
m=5.7;
kh=kf*((k1+(m-1)*kf-(m-1)*ph1*(kf-k1))/((k1+(m-1)*kf+ph1*(kf-k1))))*((k2+(m-1)*kf-(m-1)*ph2*(kf-k2))/((k2+(m-1)*kf+ph2*(kf-k2))));
muh= muf/((1-ph1)^2.5*(1-ph2)^2.5);
rhoh=rhof*(1-ph2)*((1-ph1)+ph1*(rho1/rhof))+ph2*rho2;
v1 =rhof*cpf*(1-ph2)*((1-ph1)+ph1*((rho1*cp1)/(rho2*cp2)))+ph2*(rho2*cp2);
sigh=sigf+(3*((ph1*sig1+ph2*sig2)-sigf*(ph1+ph2))/(((ph1*sig1+ph2*sig2)/(sigf*(ph1+ph2)))+2-((ph1*sig1+ph2*sig2)/sigf)+(ph1+ph2)));
alfah=kh/v1;
myLegend1 = {};
rr = [4 5 6 7]
for i =1:numel(rr)
Prf = rr(i);
Nr=0.1;
gamma=pi/3;
a=1;b=0.1;v=1;u=1;
M=3;
Nt=1;Nb=1; sc=0.6;s1=1;s2=1;
p=-0.5; L=(muf/rhof);L1=L^(p);
Tw=273+50;Ti=273+27;deltaT=Tw-Ti;
Lf=rhof*kf;
y0 = [1,0,1,0,0,1,0,1,0];
options =bvpset('stats','on','RelTol',1e-4);
m = linspace(0,4);
solinit = bvpinit(m,y0);
sol= bvp4c(@projfun,@projbc,solinit,options);
disp((sol.y(1,20)))
figure(1)
plot(sol.x,(sol.y(6,:)))
% axis([0 4 0 1])
grid on,hold on
myLegend1{i}=['Pr = ',num2str(rr(i))];
i=i+1;
end
figure(1)
legend(myLegend1)
hold on
function dy= projfun(~,y)
dy= zeros(9,1);
% alignComments
E = y(1);
dE = y(2);
F = y(3);
dF= y(4);
W = y(5);
t = y(6);
dt = y(7);
phi = y(8);
dphi = y(9);
dy(1) = dE;
dy(2) = (rhoh/muh)*((((a*u)/L1^(2)))*E^2+(1/L1)*W*dE+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*E*sin(gamma)*sin(gamma));
dy(3) = dF;
dy(4) = (rhoh/muh)*((((b*v)/L1^(2)))*F^2+(1/L1)*W*dF+((sigh/sigf)/(rhoh/rhof))*(1/L1^2)*M*F*sin(gamma)*sin(gamma));
dy(5) = -(1/L1)*(u*a*E+b*v*F);
dy(6) = dt;
dy(7) =(Prf*(rhof/muf))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2);
dy(8)= dphi;
dy(9)=(sc/L^(p+1))*W*dphi-(s1/s2)*(Nt/Nb)*(((Prf*(rhof/muf)))*(1/(Nr+(kh/kf)))*(((v1)/(rhof*cpf))*(1/L1)*W*dt-(muh/(rhof*cpf))*(L1/s1)*(1/deltaT)*(-(1/L1)*(u*a*E+b*v*F))^2));
end
end
function res= projbc(ya,yb)
res= [ya(1)-1;
ya(3)-1;
ya(5)-0;
ya(6)-1;
ya(8)-1;
yb(1);
yb(3);
yb(6);
yb(8)];
end
0 个评论
采纳的回答
Abhinaya Kennedy
2024-8-22,3:24
编辑:Walter Roberson
2024-8-22,3:40
To create a surface plot using the bvp4c solution in MATLAB, you need to organize your data into matrices that represent the x, y, and z axes. In your case, the x axis is represented by sol.x, the y axis by Prf, and the z axis by sol.y(6,:). Here's how you can modify your code to produce a surface plot:
proj()
function sol = proj
clc; clf; clear;
% Define constants and parameters
rhof = 1; kf = 0.613e5; cpf = 4179e4; muf = 1e-3 * 10; sigf = 0.05e-8;
alfaf = kf / (rhof * cpf);
ph1 = 0.01; rho1 = 5180e-3; cp1 = 670e4; k1 = 9.7e5; sig1 = 0.74e-2;
ph2 = 0.01; rho2 = 8933e-3; cp2 = 385e4; k2 = 401e5; sig2 = 5.96e-1;
m = 5.7;
kh = kf * ((k1 + (m - 1) * kf - (m - 1) * ph1 * (kf - k1)) / ((k1 + (m - 1) * kf + ph1 * (kf - k1)))) * ...
((k2 + (m - 1) * kf - (m - 1) * ph2 * (kf - k2)) / ((k2 + (m - 1) * kf + ph2 * (kf - k2))));
muh = muf / ((1 - ph1)^2.5 * (1 - ph2)^2.5);
rhoh = rhof * (1 - ph2) * ((1 - ph1) + ph1 * (rho1 / rhof)) + ph2 * rho2;
v1 = rhof * cpf * (1 - ph2) * ((1 - ph1) + ph1 * ((rho1 * cp1) / (rho2 * cp2))) + ph2 * (rho2 * cp2);
sigh = sigf + (3 * ((ph1 * sig1 + ph2 * sig2) - sigf * (ph1 + ph2)) / ...
(((ph1 * sig1 + ph2 * sig2) / (sigf * (ph1 + ph2))) + 2 - ((ph1 * sig1 + ph2 * sig2) / sigf) + (ph1 + ph2)));
alfah = kh / v1;
% Initialize parameters
rr = [4 5 6 7];
numPrf = numel(rr);
m = linspace(0, 4);
y0 = [1, 0, 1, 0, 0, 1, 0, 1, 0];
options = bvpset('stats', 'on', 'RelTol', 1e-4);
% Initialize storage for surface plot
Z = zeros(numPrf, length(m));
% Solve the BVP for each Prf
for i = 1:numPrf
Prf = rr(i);
solinit = bvpinit(m, y0);
sol = bvp4c(@projfun, @projbc, solinit, options);
Z(i, :) = sol.y(6, :); % Store the z-axis data
end
% Create surface plot
[X, Y] = meshgrid(m, rr);
figure;
surf(X, Y, Z);
xlabel('x');
ylabel('Prf');
zlabel('Solution y(6,:)');
title('Surface Plot of Solution');
grid on;
function dy = projfun(~, y)
dy = zeros(9, 1);
E = y(1);
dE = y(2);
F = y(3);
dF = y(4);
W = y(5);
t = y(6);
dt
end
end
3 个评论
更多回答(0 个)
另请参阅
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!