Could anyone help me with my matlab bvp4c program?

4 次查看（过去 30 天）
sunitha kolasani
sunitha kolasani 2024-4-30
评论： sunitha kolasani 2024-9-24
Hello Dear Sir/ Medam I have two couppled equations (temparature and Concentration) with Boundary conditions. I have solved this equations and solved analytically and got the graph. Now i have used BVP4C for the same equations and got graph is not maching with my earlier graph. Both program and graph i am sharing hear, if u help me to find where i got wrong. It will be helpful to me. Thank you.

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

采纳的回答

Torsten
Torsten 2024-4-30
Ran in:
The order of your functions in the dydy function handle is wrong.
It should work now.
%I am a Researcher, my problem is i try to solve couppled ODE s with BVP4C MATLAB code the graph obtained after running the MATLAB program, which differs from the actual graph obtained using Mathematica software. Here i am sharing the Equations and Boundary Conditions.
% temp Eqn: Theta^''=-N1*Pr*G6* (dTheta/dy)+S1*Pr*G7*Theta;
% Concen Eqn: Phi^''= -N1*Sc*G8*(dPhi/dy)+H1*Sc*G8*Phi;
% Boundary conditions: dTheta/dy=G7*B1*( Theta-1),dPhi/dy=G8*F1*( Phi-1) at y=0.
% Theta,Phi = 0 at y tends to infinity.
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.15;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
dydx = @(x,y)[y(2);
-N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
y(4);
-N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
%dydx = @(x,y)[y(2);
% y(4);
% -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
% -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
yb(1);
ya(4)- F1.*G8*(ya(3)-1);
yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T1 = bvp4c(dydx, BC, solint);
hold on
plot(T1.x,([T1.y(1,:)]),'b','linewidth',1.5)
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.25;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
% y(4);
% -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
% -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
-N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
y(4);
-N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
yb(1);
ya(4)- F1.*G8*(ya(3)-1);
yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T2 = bvp4c(dydx, BC, solint);
hold on
plot(T2.x,([T2.y(1,:)]),'r','linewidth',1.5)
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.35;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
% y(4);
% -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
% -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
-N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
y(4);
-N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
yb(1);
ya(4)- F1.*G8*(ya(3)-1);
yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T3 = bvp4c(dydx, BC, solint);
hold on
plot(T3.x,([T3.y(1,:)]),'g','linewidth',1.5)
  5 个评论
显示 3更早的评论
Torsten
Torsten 2024-9-21
Ran in:
Use more points for plotting in x-direction:
%I am a Researcher, my problem is i try to solve couppled ODE s with BVP4C MATLAB code the graph obtained after running the MATLAB program, which differs from the actual graph obtained using Mathematica software. Here i am sharing the Equations and Boundary Conditions.
% temp Eqn: Theta^''=-N1*Pr*G6* (dTheta/dy)+S1*Pr*G7*Theta;
% Concen Eqn: Phi^''= -N1*Sc*G8*(dPhi/dy)+H1*Sc*G8*Phi;
% Boundary conditions: dTheta/dy=G7*B1*( Theta-1),dPhi/dy=G8*F1*( Phi-1) at y=0.
% Theta,Phi = 0 at y tends to infinity.
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.15;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
dydx = @(x,y)[y(2);
-N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
y(4);
-N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
%dydx = @(x,y)[y(2);
% y(4);
% -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
% -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
yb(1);
ya(4)- F1.*G8*(ya(3)-1);
yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T1 = bvp4c(dydx, BC, solint);
hold on
xplot = linspace(1e-6,1.8,100);
y1plot = deval(xplot,T1);
plot(xplot,y1plot(1,:),'b','linewidth',1.5)
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.25;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
% y(4);
% -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
% -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
-N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
y(4);
-N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
yb(1);
ya(4)- F1.*G8*(ya(3)-1);
yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T2 = bvp4c(dydx, BC, solint);
hold on
xplot = linspace(1e-6,1.8,100);
y1plot = deval(xplot,T2);
plot(xplot,y1plot(1,:),'r','linewidth',1.5)
B1 = 0.5;
F1 = 0.5;
H1 = 0.2;
M1 = 1;
N1 = 0.35;
Pr = 6.2;
S1 = 1;
P= 0.02;%nano particle
Sc = 0.78;
a3 = 4420;%rows
a4 = 0.56;%(cp)s
a1 = 997.1;%rowf
a2 = 4179;%(cp)f
b = (a3.*a4)./(a1.*a2);%(row*cp)s/(row*cp)f
e6 = ((1 - P) + (P.*b));%e6
k1 = 0.613;%kf
k2 = 7.2;%ks
w = (k1./k2);%kf/ks
e7 = (((1 + 2.*P) + 2.*(1 - P).*w)./((1 - P) + (2 + P).*w));%e7
G6 = (e6./e7);
G7 = (1./e7);
G8 = (1./(1-P));
%dydx = @(x,y)[y(2);
% y(4);
% -N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
% -N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
dydx = @(x,y)[y(2);
-N1.*Pr.*G6.*y(2)+S1.*Pr.*G7.*y(1);%TEMP
y(4);
-N1.*Sc.*G8.*y(4)+H1.*Sc.*G8.*y(3)];%concentration
BC = @(ya,yb)[ya(2)-G7.*B1.*(ya(1)-1);
yb(1);
ya(4)- F1.*G8*(ya(3)-1);
yb(3)];
yinit = [0.1;0.1;0.1;0.1];
solint = bvpinit(linspace(1e-6,1.8,15),yinit);
T3 = bvp4c(dydx, BC, solint);
hold on
xplot = linspace(1e-6,1.8,100);
y1plot = deval(xplot,T3);
plot(xplot,y1plot(1,:),'g','linewidth',1.5)
%plot(T3.x,([T3.y(1,:)]),'g','linewidth',1.5)

请先登录，再进行评论。

更多回答（0 个）

类别

在 Help Center 和 File Exchange 中查找有关 Solver Outputs and Iterative Display 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by