清空筛选条件
清空筛选条件

Please help me to run surf plot with bvp4c.

34 次查看（过去 30 天）
T K
T K 2023-2-7
评论： Walter Roberson 2024-8-22，21:34
Ran in:
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()
rr = 1x4
4 5 6 7
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
The solution was obtained on a mesh of 227 points. The maximum residual is 8.171e-05. There were 7965 calls to the ODE function. There were 122 calls to the BC function. 0.2422 The solution was obtained on a mesh of 179 points. The maximum residual is 2.747e-05. There were 7696 calls to the ODE function. There were 149 calls to the BC function. 0.1666 The solution was obtained on a mesh of 188 points. The maximum residual is 2.582e-05. There were 7831 calls to the ODE function. There were 149 calls to the BC function. 0.1666 The solution was obtained on a mesh of 198 points. The maximum residual is 3.095e-05. There were 7981 calls to the ODE function. There were 149 calls to the BC function. 0.1935
ans = struct with fields:
solver: 'bvp4c' x: [0 0.0202 0.0404 0.0808 0.1212 0.1818 0.2424 0.3030 0.3636 0.4242 0.4848 0.5253 0.5657 0.6061 0.6465 0.7071 0.7677 0.8081 0.8485 0.8889 0.9293 ... ] (1x198 double) y: [9x198 double] yp: [9x198 double] stats: [1x1 struct]
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

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

采纳的回答

Abhinaya Kennedy
Abhinaya Kennedy 2024-8-22，3:24
编辑：Walter Roberson 2024-8-22，3:40
Ran in:
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:
  1. Initialize Data Storage: Store the results for each value of Prf in a matrix.
  2. Loop Over Prf: Solve the boundary value problem for each value of Prf.
  3. Create the Surface Plot: Use surf to plot the data.
proj()
Unrecognized function or variable 'dt'.

Error in solution>proj/projfun (line 56)
dt

Error in bvparguments (line 96)
testODE = ode(x1,y1,odeExtras{:});

Error in bvp4c (line 119)
bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in solution>proj (line 34)
sol = bvp4c(@projfun, @projbc, solinit, options);
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 个评论
显示 1更早的评论
T K
T K 2024-8-22，12:55
编辑：Walter Roberson 2024-8-22，21:17
Ran in:
proj()
rr = 1x4
4 5 6 7
<mw-icon class=""></mw-icon>
<mw-icon class=""></mw-icon>
The solution was obtained on a mesh of 227 points. The maximum residual is 8.171e-05. There were 7965 calls to the ODE function. There were 122 calls to the BC function. 0.2422 The solution was obtained on a mesh of 179 points. The maximum residual is 2.747e-05. There were 7696 calls to the ODE function. There were 149 calls to the BC function. 0.1666 The solution was obtained on a mesh of 188 points. The maximum residual is 2.582e-05. There were 7831 calls to the ODE function. There were 149 calls to the BC function. 0.1666 The solution was obtained on a mesh of 198 points. The maximum residual is 3.095e-05. There were 7981 calls to the ODE function. There were 149 calls to the BC function. 0.1935
ans = struct with fields:
solver: 'bvp4c' x: [0 0.0202 0.0404 0.0808 0.1212 0.1818 0.2424 0.3030 0.3636 0.4242 0.4848 0.5253 0.5657 0.6061 0.6465 0.7071 0.7677 0.8081 0.8485 0.8889 0.9293 ... ] (1x198 double) y: [9x198 double] yp: [9x198 double] stats: [1x1 struct]
%full code is showen below...
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
Walter Roberson
Walter Roberson 2024-8-22，21:34
Okay... but you do not produce a surface plot?

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 General Physics 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by