PDE solve with bvp4c
4 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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.
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1623693/image.jpeg)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/1623698/image.jpeg)
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);
0 Kommentare
Antworten (1)
Torsten
am 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:
6 Kommentare
Torsten
am 23 Feb. 2024
Bearbeitet: Torsten
am 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
Siehe auch
Kategorien
Mehr zu Boundary Value Problems finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!