What changes are needed to run the code
Afficher commentaires plus anciens
Psi = pi/2; L1 = 0.1; L2 = 0.1; L = 0.1; Pr = 1; M = 5;
Ec = 0.5; Nb = 0.5; Nt = 0.1; fw = 0.5; Q = 0.1; Le = 2; Kc = 1;
A = linspace(-1,1,101);
for L = [1 2 3]/10
BC = @(ya,yb)[ ya(1)-fw; ya(2)+1; ya([4;7])-1; yb([2;4;6]) ];
ODE = @(x,y)[ y(2); y(3); ( A.*(y(2)+(1/2)*x*y(3)) - y(1)*y(3) + y(2)*( y(2)+M^2*(sin(Psi))^2) - L1*y(4) - L2*y(6))/(1+(1/L));
y(5); -Pr*( y(1)*y(5) - (1/2)*A.*x*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2 + M^2*Ec*y(2)^2+Q*y(4) );
y(7); -Le*Pr*( y(1)*y(7)-(1/2)*A.*x*y(7)-Kc*y(6) ) - (Nt/Nb)*(-Pr*( y(1)*y(5) - (1/2)*A.*x*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2 + M^2*Ec*y(2)^2+Q*y(4) )) ];
xa = 0; xb = 6; xn = linspace(xa,xb,101); x = xn; solinit = bvpinit(x,[0 1 0 1 0 1 0]); sol = bvp5c(ODE,BC,solinit);
S = deval(sol,x); f0 = deval(sol,xa);
figure(1);plot(A,S(2,:),'-','LineWidth',2);xlabel('\bf\eta','color','b'); ylabel('\bff^\prime(\eta)','color','b');hold on,
end
%% ERROR is:
Dimensions of matrices being concatenated are not consistent.
%% Help needed to run
Réponses (1)
ODE = @(x,y) [ ...
y(2); ...
y(3); ...
(A.*(y(2)+0.5*x*y(3)) - y(1)*y(3) + y(2)*( y(2)+M^2*(sin(Psi))^2) - L1*y(4) - L2*y(6))/(1+(1/L)); ...
y(5); ...
-Pr*(y(1)*y(5) - 0.5*A.*x*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2 + M^2*Ec*y(2)^2+Q*y(4)); ...
y(7); ...
-Le*Pr*(y(1)*y(7) - 0.5*A.*x*y(7)-Kc*y(6)) - (Nt/Nb)*(-Pr*(y(1)*y(5) - 0.5*A.*x*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2 + M^2*Ec*y(2)^2+Q*y(4)))];
Yes, the 3rd, 5th and 7th element have the size 1x101, while the others are scalars. This comes from the vector A. I cannot guess, what you want to do instead.
7 commentaires
MINATI PATRA
le 16 Nov 2021
Then you have to run the integration for all values of A. You cannot find all different trajectories with one integration.
Psi = pi/2;
L1 = 0.1;
L2 = 0.1;
L = 0.1;
Pr = 1;
M = 5;
Ec = 0.5;
Nb = 0.5;
Nt = 0.1;
fw = 0.5;
Q = 0.1;
Le = 2;
Kc = 1;
A = linspace(-1, 1, 10);
xa = 0;
xb = 6;
xn = linspace(xa,xb,101);
x = xn;
figure;
axesH = axes('NextPlot', 'add');
for L = [1 2 3]/10
for iA = A
BC = @(ya,yb)[ ya(1)-fw; ya(2)+1; ya([4;7])-1; yb([2;4;6]) ];
ODE = @(x,y) [ ...
y(2); ...
y(3); ...
(iA.*(y(2)+0.5*x*y(3)) - y(1)*y(3) + y(2)*( y(2)+M^2*(sin(Psi))^2) - L1*y(4) - L2*y(6))/(1+(1/L)); ...
y(5); ...
-Pr*(y(1)*y(5) - 0.5*iA.*x*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2 + M^2*Ec*y(2)^2+Q*y(4)); ...
y(7); ...
-Le*Pr*(y(1)*y(7) - 0.5*iA.*x*y(7)-Kc*y(6)) - (Nt/Nb)*(-Pr*(y(1)*y(5) - 0.5*iA.*x*y(5) + Nb*y(5)*y(7) + Nt*y(5)^2 + M^2*Ec*y(2)^2+Q*y(4)))];
solinit = bvpinit(x, [0 1 0 1 0 1 0]);
sol = bvp5c(ODE, BC, solinit);
S = deval(sol, x);
f0 = deval(sol, xa);
plot(axesH, iA, S(2,:), '.', 'LineWidth', 2);
end
end
xlabel('\bf\eta','color','b');
ylabel('\bff^\prime(\eta)','color','b');
I assume you want to modify the oputput. Maybe collect the different output of S(:,2) in an array at first and draw it as plot(A, collectedS)?
MINATI PATRA
le 21 Nov 2021
Jan
le 22 Nov 2021
Please use the tools in the forum to format your code. This improves the readability.
Why do you expect 3 curves only?
MINATI PATRA
le 22 Nov 2021
MINATI PATRA
le 6 Déc 2021
Jan
le 6 Déc 2021
I'm not sure, which problem is still open. The code of my answer is running, and only the diagram looks a little bit strange. I do not know, what you want to get instead of this solution. So please explain again, what is still unsolved.
Catégories
En savoir plus sur Programming dans Centre d'aide et File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
