# PDE solve with bvp4c

15 views (last 30 days)
mallela ankamma rao on 21 Feb 2024
Edited: Torsten on 23 Feb 2024
Good evening Torsten sir
I am trying to solve the system of eight partial differential equations using bvp4c solver. I got the following error. I request you please help me to get correct graph. I also attached my equations pics sir.
MY CODE :
syms y1(t) y2(t) y3(t) y4(t) y5(t) y6(t) y7(t) y8(t)
w=1;a=1;m=2;k1=1;k2=1;k3=1;k4=1;k=1; ps=1;po=1;r1=0;r2=2;
dy1=diff(y1,t);
dy2=diff(y2,t);
dy3=diff(y3,t);
dy4=diff(y4,t);
dy5=diff(y5,t);
dy6=diff(y6,t);
dy7=diff(y7,t);
dy8=diff(y8,t);
eq1=diff(y1,2)==(m^2+(1/k1))*y1-ps*r1;
eq2=diff(y2,2)==(m^2+(1/k2))*y2+a*ps*r2;
eq3=diff(y3,2)==-k1*(diff(y1,1))^2;
eq4=diff(y4,2)==-k2*(diff(y2,1))^2;
eq5=diff(y5,2)==(m^2+(1/k1)+1i*w*r1)*y5-po*r1;
eq6=diff(y6,2)==(m^2+(1/k1)+1i*w*r2)*y6-a*po*r2;
eq7=diff(y7,2)==k*y7-k3*diff(y1,1)*diff(y3,1);
eq8=diff(y8,2)==k*y8-k4*diff(y2,1)*diff(y4,1);
vars = [y1(t); y2(t); y3(t); y4(t);y5(t); y6(t); y7(t); y8(t)];
V = odeToVectorField([eq1,eq2,eq3, eq4, eq5,eq6,eq7,eq8]);
M = matlabFunction(V,'vars', {'t','Y'});
t = linspace(0,10,100);
init = bvpinit(t,[0,0,1,1,0,0,0,0]);
%Run ODE solver
options = bvpset('RelTol',1e-8,'AbsTol',1e-12);
sol = bvp4c(M,@bcfun,init,options);
% plot(sol.y, sol.t)
function res = bcfun(yl,y2,y3,y4,y5,y6,y7,y8,dyl,dy2,dy3,dy4,dy5,dy6,dy7,dy8)
%res=[yl(1)-1; y2(-1); y1(0)-y2(0); u1*dyl(0)-u2*dy2(0); y5(1)-1; y6(-1);y5(0)-y6(0);u1*dy5(0)-u2*dy6(0);y3(-1);y4(1)-1;y3(0)-y4(0);k1*dy3(0)-k2*dy4(0);y7(-1);y8(1);y7(0)-y8(0);k1*dy7(0)-k2*dy8(0)];
res = zeros(16,1);
res(1) = yl(1)-1;
res(2) = y2(-1);
res(3) = y1(0)-y2(0);
res(4) = u1*dyl(0)-u2*dy2(0);
res(5) = y5(1)-1;
res(6) = y6(-1);
res(7) = y5(0)-y6(0);
res(8) = u1*dy5(0)-u2*dy6(0);
res(9) = y3(-1);
res(10) = y4(1)-1;
res(11) = y3(0)-y4(0);
res(12) = k1*dy3(0)-k2*dy4(0);
res(13) = y7(-1);
res(14) = y8(1);
res(15) = y7(0)-y8(0);
res(16) = k1*dy7(0)-k2*dy8(0);
end
plot(sol.t,sol.y(2,:),'-')
ERROR:
Index exceeds the number of array elements (8).
Error in
symengine>@(t,Y)[Y(2);Y(1).*5.0+2.0;Y(4);Y(3).*5.0;Y(6);-Y(4).^2;Y(8);-Y(2).^2;Y(10);Y(9).*5.0;Y(12);Y(11).*(5.0+2.0i)-2.0;Y(14);Y(13)-Y(4).*Y(6);Y(16);Y(15)-Y(2).*Y(8)]
Error in bvparguments (line 105)
testODE = ode(x1,y1,odeExtras{:});
Error in bvp4c (line 128)
bvparguments(solver_name,ode,bc,solinit,options,varargin);
Error in trial (line 34)
sol = bvp4c(M,@bcfun,init,options);

Torsten on 21 Feb 2024
Take a look here on how boundary value problems with two regions (in your case -1 to 0 and 0 to 1) have to be set up for bvp4c:
Torsten on 23 Feb 2024
Edited: Torsten on 23 Feb 2024
I think you should be able to add the parameters needed.
I'm not sure whether the part
dydx(1) = usx;
dydx(2) = (M^2 + 1/K2)*us - alpha*Ps*R2;
dydx(3) = uox;
dydx(4) = (M^2 + 1/K1 + 1i*omega*R2)*uo - alpha*Po*R2; % 1/K2 ?
belongs to region 2 (because of K2 and R2) and the part
dydx(1) = usx;
dydx(2) = (M^2 + 1/K1)*us - Ps*R1;
dydx(3) = uox;
dydx(4) = (M^2 + 1/K1 + 1i*omega*R1)*uo - Po*R1;
belongs to region 1 (because of K1 and R1)
But you prescribed boundary conditions for u21 and u22 at y=-1 and for u11 and u12 at y=1. Thus this looks a little contradictory.
xc = 0;
npts = 20;
xmesh1 = linspace(-1,xc,npts);
xmesh2 = linspace(xc,1,npts);
xmesh = [xmesh1,xmesh2];
yinit = [0 1 0 1 0 1 0 1];
init = bvpinit(xmesh,yinit);
sol = bvp4c(@f, @bc, init);
function dydx = f(x,y,region) % equations being solved
us = y(1);
usx = y(2);
uo = y(3);
u0x = y(4);
thetas = y(5);
thetasx = y(6);
thetao = y(7);
thetaox = y(8);
dydx = zeros(8,1);
switch region
case 1 % x in [-1 0]
dydx(1) = usx;
dydx(2) = (M^2 + 1/K2)*us - alpha*Ps*R2;
dydx(3) = uox;
dydx(4) = (M^2 + 1/K1 + 1i*omega*R2)*uo - alpha*Po*R2; % 1/K2 ?
dydx(5) = thetasx;
dydx(6) = k1*usx^2;
dydx(7) = thetaox;
dydx(8) = k*thetao - k3*thetasx*thetaox;
case 2 % x in [0 1]
dydx(1) = usx;
dydx(2) = (M^2 + 1/K1)*us - Ps*R1;
dydx(3) = uox;
dydx(4) = (M^2 + 1/K1 + 1i*omega*R1)*uo - Po*R1;
dydx(5) = thetasx;
dydx(6) = k2*usx^2;
dydx(7) = thetaox;
dydx(8) = k*thetao - k4*thetasx*thetaox;
end
end
function res = bc(YL,YR)
res = [YL(1,1)
YL(3,1)
YL(5,1)
YL(7,1)
YR(1,2)-1
YR(3,2)-1
YR(5,2)-1
YR(7,2)
YR(1,1)-YL(1,2)
mu1*YR(2,1)-mu2*YL(2,2)
YR(3,1)-YL(3,2)
mu1*YR(4,1)-mu2*YL(4,2)
YR(5,1)-YL(5,2)
k1*YR(6,1)-k2*YL(6,2)
YR(7,1)-YL(7,2)
k1*YR(8,1)-k2*YL(8,2)];
end
mallela ankamma rao on 23 Feb 2024
Thank you verymuch for spending your valuable time Torsten sir.
I am truly appreciative of your assistance, sir.