Equations and Boundary conditions are Unequal
2 次查看(过去 30 天)
显示 更早的评论
%%% THe present code is of the attached (Model#02) pdf, need modification to run.
%%% REFERENCE: This type of work was done in DOI: 10.1002/num.22672 (More in size, so cant be uploaded)
main()
function main
We = 2; n = 0.5; M = 1; Pr = 7; Rd = 0.5; Tw = 1.5; Nt = 0.3; Nb = 0.5; Ec = 0.3; Q = 0.5; Le = 2; K = 0.1; D1 = 0.1;
m = 0.5; E = 0.1; Lb = 0.5; Pe = 2; D2 = 1; S1 = 0.2; S2 = 0.2; Om = 0.5;
xa = 0; xb = 5; solinit = bvpinit(linspace(xa,xb,100),[0 1 0 1 0 1 0 1 0 1 0 1]); sol = bvp5c(@ode,@bc,solinit); x = linspace(xa,xb,100); S = deval(sol,x);
function BC = bc(ya,yb)
BC = [ya([1,4]); ya(2) - S1; ya([6,8,10]) - 1; yb([6,8,10]); yb(2) - S2; yb(4) - Om];
end
function EQ = ode(x,y)
Af = (1+(We*y(3))^2)^((n-1)/2) + (n-1)*(We*y(3))^2*(1+(We*y(3))^2)^((n-3)/2); Ag = (1+(We*y(5))^2)^((n-1)/2) + (n-1)*(We*y(5))^2*(1+(We*y(5))^2)^((n-3)/2);
At = 4*Rd*(Tw-1)*y(7)^2*(1+(Tw-1)*y(6))^2; X = - E/(1+D1*y(6));
EQ = [ -2*y(2);
y(3); (M*y(2) + y(2)^2 - y(4)^2 + y(1)*y(3))/Af;
y(5); (M*y(4) + 2*y(2)*y(4) + y(1)*y(5))/Ag;
y(7); (Pr/At)*( y(1)*y(7) - Nt*y(7)^2 - Nb*y(7)*y(9) - Ec*(y(3)^2 + y(5)^2) - M*Ec*(y(2)^2 + y(4)^2) - Q*y(6) );
y(9); Pr*Le*( y(1)*y(9) + K*(1+D1*y(6))^m*y(8)*exp(X) - (Nt/Nb)*((Pr/At)*( y(1)*y(7) - Nt*y(7)^2 - Nb*y(7)*y(9) - Ec*(y(3)^2 + y(5)^2) - M*Ec*(y(2)^2 + y(4)^2) - Q*y(6) )) );
y(11); Lb*y(1)*y(11)+ Pe*( y(9)*y(11) + (D2 + y(10)) )*Pr*Le*( y(1)*y(9) + K*(1+D1*(y(6))^m)*y(8)*exp(X) - (Nt/Nb)*((Pr/At)*( y(1)*y(7) - Nt*y(7)^2 - Nb*y(7)*y(9) - Ec*(y(3)^2 + y(5)^2) - M*Ec*(y(2)^2 + y(4)^2) - Q*y(6) ) ))
];
figure(10),plot(x,S(2,:),'-b','LineWidth',1.5),hold on
end
figure(2),plot(x,S(2,:));hold on
end
0 个评论
采纳的回答
Walter Roberson
2023-8-14
移动:Walter Roberson
2023-8-14
S is the output of the deval() and so is not available until after the bvp5c has been run. But the ode function wants to plot S, which requires S be defined but it is not defined until after the ode is finished.
Your plotting should be moved to before the bc function definition.
Meanwhile, your ode function needs to return a column vector of length 12, same length as your initial condition; at the moment it is only length 11.
main()
function main
We = 2; n = 0.5; M = 1; Pr = 7; Rd = 0.5; Tw = 1.5; Nt = 0.3; Nb = 0.5; Ec = 0.3; Q = 0.5; Le = 2; K = 0.1; D1 = 0.1;
m = 0.5; E = 0.1; Lb = 0.5; Pe = 2; D2 = 1; S1 = 0.2; S2 = 0.2; Om = 0.5;
xa = 0; xb = 5;
solinit = bvpinit(linspace(xa,xb,100),[0 1 0 1 0 1 0 1 0 1 0 1]);
sol = bvp5c(@ode,@bc,solinit);
x = linspace(xa,xb,100);
S = deval(sol,x);
figure(10),plot(x,S(2,:),'-b','LineWidth',1.5),hold on
function BC = bc(ya,yb)
BC = [ya([1,4]); ya(2) - S1; ya([6,8,10]) - 1; yb([6,8,10]); yb(2) - S2; yb(4) - Om];
end
function EQ = ode(x,y)
Af = (1+(We*y(3))^2)^((n-1)/2) + (n-1)*(We*y(3))^2*(1+(We*y(3))^2)^((n-3)/2); Ag = (1+(We*y(5))^2)^((n-1)/2) + (n-1)*(We*y(5))^2*(1+(We*y(5))^2)^((n-3)/2);
At = 4*Rd*(Tw-1)*y(7)^2*(1+(Tw-1)*y(6))^2; X = - E/(1+D1*y(6));
EQ = [ -2*y(2);
y(3); (M*y(2) + y(2)^2 - y(4)^2 + y(1)*y(3))/Af;
y(5); (M*y(4) + 2*y(2)*y(4) + y(1)*y(5))/Ag;
y(7); (Pr/At)*( y(1)*y(7) - Nt*y(7)^2 - Nb*y(7)*y(9) - Ec*(y(3)^2 + y(5)^2) - M*Ec*(y(2)^2 + y(4)^2) - Q*y(6) );
y(9); Pr*Le*( y(1)*y(9) + K*(1+D1*y(6))^m*y(8)*exp(X) - (Nt/Nb)*((Pr/At)*( y(1)*y(7) - Nt*y(7)^2 - Nb*y(7)*y(9) - Ec*(y(3)^2 + y(5)^2) - M*Ec*(y(2)^2 + y(4)^2) - Q*y(6) )) );
y(11); Lb*y(1)*y(11)+ Pe*( y(9)*y(11) + (D2 + y(10)) )*Pr*Le*( y(1)*y(9) + K*(1+D1*(y(6))^m)*y(8)*exp(X) - (Nt/Nb)*((Pr/At)*( y(1)*y(7) - Nt*y(7)^2 - Nb*y(7)*y(9) - Ec*(y(3)^2 + y(5)^2) - M*Ec*(y(2)^2 + y(4)^2) - Q*y(6) ) ))
];
size(EQ)
end
figure(2),plot(x,S(2,:));hold on
end
4 个评论
Torsten
2023-8-15
You generate an artificial degree of freedom by this differentiation. Depending on the second boundary condition you impose you may or may not reproduce the solution of the original problem (first-order ODE with only one boundary condition).
Test it for a simple example.
更多回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Ordinary Differential Equations 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!