Please help me to run surf plot with bvp4c.

Question

T K 2023-2-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1907905-please-help-me-to-run-surf-plot-with-bvp4c

评论： T K 2024-9-1

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

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

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

请先登录，再进行评论。

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

Answer 1

Abhinaya Kennedy 2024-8-22

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1907905-please-help-me-to-run-surf-plot-with-bvp4c#answer_1502939

编辑：Walter Roberson 2024-8-22

在 MATLAB Online 中打开

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:

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.

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

7 个评论
显示 5更早的评论隐藏 5更早的评论

T K 2024-8-22

编辑：Walter Roberson 2024-8-22

在 MATLAB Online 中打开

proj()

rr = 1x4

4 5 6 7

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

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

T K 2024-9-1

在 MATLAB Online 中打开

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);
    a=1;a=1;b=0.1;v=1;u=1;
    M=3;    gamma=pi/3;
Nr=0.1;
     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);
    
    % 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);
        % 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

Torsten 2024-9-1

在 MATLAB Online 中打开

It may happen that "bvp4c" changes the number of discretizations points so that it won't return the solution in 100 points as at the start, but in more or less.

Use

sol = bvp4c(@projfun, @projbc, solinit, options);
Z(i,:) = deval(sol,m,6)

to interpolate the solution to the values m of your choice.

Take care not to come into conflict with the parameter "m" used here:

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))));

T K 2024-9-1

the m in the relation kh replaced by s s=5.7, Kh=.....

请先登录，再进行评论。

Please help me to run surf plot with bvp4c.

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论
显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

Please help me to run surf plot with bvp4c.

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

采纳的回答

7 个评论 显示 5更早的评论隐藏 5更早的评论

更多回答（0 个）

另请参阅

类别

标签

Community Treasure Hunt

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

7 个评论
显示 5更早的评论隐藏 5更早的评论