Filter löschen
Filter löschen

Getting error while trying to solve boundary condition equations.

7 Ansichten (letzte 30 Tage)
krishna juttiga
krishna juttiga am 19 Apr. 2024
Beantwortet: Torsten am 19 Apr. 2024
I'm currently working on solving a boundary value problem (BVP) represented by a system of ordinary differential equations in MATLAB. However, I've encountered some errors during the process, and I'm seeking assistance in resolving them. Below is the code snippet
close all; clear all; clc;
K1=0; eta1=0.2; eta2=0.7; lambda1=0.5; Fr=0.5; Rd=0.5; Sc=0.8; Pr=1; omega=0.5; Nr=0.5; B1=0.5; Nt=0.4; beta1=0.3; m1=0.4; Pe=0.3; E2=0.2; Lb=0.4; Rb=0.6; alpha=0.3; r1=0.3; Kp=0.3;
w=1; Nb=0.5; B2=0.5;
% fprintf('Working\n');
for M1 = 0:0.1:1
V1(w)= M1;
j=1;
% fprintf('Working\n');
% global khp;
% khp =@ (x,y) [y(2);
% y(3);
% (-y(1)*y(3)+(1+Fr)*(y(2)^2)+Kp*y(2)+M1*y(2)-lambda1*(y(5)-Nr*y(8)-Rb*y(11)))/(1+K1-eta1*K1*(y(3)^2));
% y(6);
% -(Pr*y(6)*y(1)+Pr*Nb*y(6)*y(9)+Pr*Nt*((y(6))^2)+alpha*Pr*y(5))/(1+Rd*4/3);
% y(9);
% (Sc*beta1*(((1+y(5)*r1))^M1)*y(8)*exp(-E2/(1+r1*y(5)))-Sc*y(1)*y(9)-(Nt/Nb)*y(7));
% y(12);
% (Pe*(y(11)+omega)*y(10)-y(12)*(Lb*y(1)-Pe*y(9)))];
global khp;
khp =@ (x,y) [y(2);
y(3);
((-y(1)*y(3))+((1+Fr)*(y(2)^2))+(Kp*y(2))+(M1*y(2))-(lambda1*(y(4)-Nr*y(6)-Rb*y(8))))/((1+K1)-(eta1*K1*(y(3)^2)));
y(5);
(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)));
y(7);
(Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)))));
y(9);
(Pe*(y(8)+omega)*((Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3))))))+((y(9))*((Pe*y(6))-(Lb*y(1)))))];
% fprintf('Working\n');
init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0]);
% fprintf('Working\n');
bc_output = bc_bvp(sol1);
sol1=bvp4c(khp,bc_output,init);
% fprintf('Working\n');
x=linspace(0,9,100);
y1=deval(sol1,x);
Cf(w,j)= ((1+K1)*y1(3,1)-((K1*eta1)/3)*y1(3,1)^(3));
Nu(w,j)=-(1+Rd*4/3)*y1(5,1);
Sh(w,j)=-(y1(7,1));
Nn(w,j)=-(y1(9,1));
w=w+1;
end
%
M1 = 0:0.1:1;
plot(M1,Cf)
function res = bc_bvp(yl,yr) % boundary conditions
eta1=0.2; eta2=0.7;K1=0;B1=0.5; B2=0.5;
%global khp;
%init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0, 0,0,0]);
bc_output = bc_bvp(sol1);
sol1=bvp4c(khp,bc_output,init);
%y1=deval(sol1,x);
% fprintf("%f\n",length(yr));
%fprintf("%f\n",yr(11));
res=[yl(1); yl(2)-1-(eta2)*((1+K1)*yl(3)-(eta1)*K1*(yl(3)^3)); yl(5)+B1*(1-yl(4)); yl(7)+B2*(1-yl(6)); yl(9)-1; yr(2); yr(4); yr(6); yr(9);];
% fprintf("%f\n",length(res));
end
looking for some guidance or suggestions .

Antworten (2)

Pavan Sahith
Pavan Sahith am 19 Apr. 2024
It appears that you are encountering an issue while trying to solve a boundary value problem (BVP) using MATLAB. The error you're experiencing is due to the use of the variable ‘sol1’ before it has been properly defined. To address this problem, you need to ensure that ‘sol1’ is defined by calling “bvp4c” before attempting to access or use it.
Issue might be solved by passing “bc_bvp” as a function handle to “bvp4c” instead of using the 'sol1' before defining.
In MATLAB, the “bvp4c” function is designed to solve BVPs and requires three key arguments:
  1. The function that defines the differential equations (‘khp’ in your case).
  2. The function that specifies the boundary conditions (‘bc_bvp’ in your case).
  3. An initial guess for the solution structure (‘init’ in your case).
You can pass the boundary condition function “bc_bvp” as a function handle using the @ symbol. The correct way to call ‘bvp4c” and define ‘sol1’ would be as follows:
sol1 = bvp4c(khp, @bc_bvp, init);
For additional guidance and examples on how to use “bvp4c”, including how to pass parameters to the function, you can refer to the MathWorks documentation.
Hope this will help you in solving the error.

Torsten
Torsten am 19 Apr. 2024
For M1 = 0-0.4, bvp4c cannot find a solution. So I started with M1 = 0.5.
close all; clear all; clc;
K1=0; eta1=0.2; eta2=0.7; lambda1=0.5; Fr=0.5; Rd=0.5; Sc=0.8; Pr=1; omega=0.5; Nr=0.5; B1=0.5; Nt=0.4; beta1=0.3; m1=0.4; Pe=0.3; E2=0.2; Lb=0.4; Rb=0.6; alpha=0.3; r1=0.3; Kp=0.3;
w=1; Nb=0.5; B2=0.5;
% fprintf('Working\n');
for M1 = 0.5:0.1:1
V1(w)= M1;
%j=1;
% fprintf('Working\n');
% global khp;
% khp =@ (x,y) [y(2);
% y(3);
% (-y(1)*y(3)+(1+Fr)*(y(2)^2)+Kp*y(2)+M1*y(2)-lambda1*(y(5)-Nr*y(8)-Rb*y(11)))/(1+K1-eta1*K1*(y(3)^2));
% y(6);
% -(Pr*y(6)*y(1)+Pr*Nb*y(6)*y(9)+Pr*Nt*((y(6))^2)+alpha*Pr*y(5))/(1+Rd*4/3);
% y(9);
% (Sc*beta1*(((1+y(5)*r1))^M1)*y(8)*exp(-E2/(1+r1*y(5)))-Sc*y(1)*y(9)-(Nt/Nb)*y(7));
% y(12);
% (Pe*(y(11)+omega)*y(10)-y(12)*(Lb*y(1)-Pe*y(9)))];
%global khp;
khp =@ (x,y) [y(2);
y(3);
((-y(1)*y(3))+((1+Fr)*(y(2)^2))+(Kp*y(2))+(M1*y(2))-(lambda1*(y(4)-Nr*y(6)-Rb*y(8))))/((1+K1)-(eta1*K1*(y(3)^2)));
y(5);
(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)));
y(7);
(Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3)))));
y(9);
(Pe*(y(8)+omega)*((Sc*beta1*(((1+(y(4)*r1)))^m1)*y(6)*exp(-E2/(1+(r1*y(4))))-(Sc*y(1)*y(7))-((Nt/Nb)*(-((Pr*y(5)*y(1))+(Pr*Nb*y(5)*y(7))+(Pr*Nt*((y(5))^2))+(alpha*Pr*y(4)))/(1+(Rd*4/3))))))+((y(9))*((Pe*y(6))-(Lb*y(1)))))];
% fprintf('Working\n');
init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0]);
% fprintf('Working\n');
%bc_output = bc_bvp(sol1);
sol1=bvp4c(khp,@bc_bvp,init);
% fprintf('Working\n');
x=linspace(0,9,100);
y1=deval(sol1,x);
Cf(w)= ((1+K1)*y1(3,1)-((K1*eta1)/3)*y1(3,1)^(3));
Nu(w)=-(1+Rd*4/3)*y1(5,1);
Sh(w)=-(y1(7,1));
Nn(w)=-(y1(9,1));
w=w+1;
end
%
M1 = 0.5:0.1:1;
plot(M1,Cf)
function res = bc_bvp(yl,yr) % boundary conditions
eta1=0.2; eta2=0.7;K1=0;B1=0.5; B2=0.5;
%global khp;
%init=bvpinit(linspace(0,9,100),[1/2, 0, 0, 1, 0, 0, 1, 0, 0, 0,0,0]);
%bc_output = bc_bvp(sol1);
%sol1=bvp4c(khp,bc_output,init);
%y1=deval(sol1,x);
% fprintf("%f\n",length(yr));
%fprintf("%f\n",yr(11));
res=[yl(1)
yl(2)-1-(eta2)*((1+K1)*yl(3)-(eta1)*K1*(yl(3)^3))
yl(5)+B1*(1-yl(4))
yl(7)+B2*(1-yl(6))
yl(9)-1
yr(2)
yr(4)
yr(6)
yr(9)];
% fprintf("%f\n",length(res));
end

Produkte


Version

R2023a

Community Treasure Hunt

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

Start Hunting!

Translated by