PDE solve with bvp4c

Question

mallela ankamma rao 2024-2-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2085018-pde-solve-with-bvp4c

编辑： Torsten 2024-2-23

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);

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Torsten 2024-2-21

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/2085018-pde-solve-with-bvp4c#answer_1413948

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:

https://uk.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html

6 个评论
显示 4更早的评论隐藏 4更早的评论

mallela ankamma rao 2024-2-23

Good evening sir

Sir I can not understand how to divides two regions and how get graphs. I humbly request you please give your valuable time to correct the code sir.

My Code:

function bvp4c

sol1 = bvpinit(linspace(-1,0,25),[1 1 0 1 1 1 1 1 1 1 1 1 1 1 1 1]);

options = bvpset('RelTol',1e-8,'AbsTol',1e-12);

sol = bvp4c(@odefun, @bcfun, sol1,options);

x = sol.x;

y = sol.y;

% Plotting of the velocity

figure (1)

plot(x, y(2, :) ,'linewidth', 2)

figure (2)

plot(x, y(4, :) ,'linewidth', 2)

% Residual of the boundary conditions

function residual = bcfun(ya, yb, yc)

residual=[yc(1)-1; ya(3); yb(1)-yb(3); u1*yb(2)-u2*yb(3); yc(9)-1; ya(11); yb(9)-yb(11); u1*dyb(10)-u2*dyb(12); ya(5); yc(7)-1; yb(5)-yb(7);...

k1*yb(6)-k2*yb(8);ya(13);yc(15);yb(13)-yb(15); k1*dyb(14)-k2*dyb(16)];

end

%System of First Order ODEs

function dy = odefun(t,y)

w=1;a=1;m=2;k1=1;k2=1;k3=1;k4=1;k=1; ps=1;po=1;r1=0;r2=2;

% y1=u11; y2=dy1; y3=u21; y4=dy3; y5=theta11; y6=dy5; y7=theta21; y8=dy7;

% y9=u12; y10=dy9; y11=u22; y12=dy11; y13=theta12; y14=dy13; y15=theta22; y16=dy15;

dy2=(m^2+(1/k1))*y(1)-ps*r1;

dy4=(m^2+(1/k2))*y(3)+a*ps*r2;

dy6=-k1*(y(2))^2;

dy8=-k2*(y(4))^2;

dy10=(m^2+(1/k1)+1i*w*r1)*y(9)-po*r1;

dy12=(m^2+(1/k1)+1i*w*r2)*y(11)-a*po*r2;

dy14=k*y(13)-k3*y(2)*y(10);

dy16=k*y(15)-k4*y(4)*y(12);

dy = [y(2);dy2;y(4);dy4;y(6);dy6;y(8);dy8;y(10);dy10;y(12);dy12;y(14);dy14;y(16);dy16];

end

error

Not enough input arguments.

Error in trial11/bcfun (line 91)

residual=[yc(1)-1; ya(3); yb(1)-yb(3); u1*yb(2)-u2*yb(3); yc(9)-1; ya(11); yb(9)-yb(11); u1*dyb(10)-u2*dyb(12); ya(5); yc(7)-1;

yb(5)-yb(7);...

Error in bvparguments (line 106)

testBC = bc(ya,yb,bcExtras{:});

Error in bvp4c (line 128)

bvparguments(solver_name,ode,bc,solinit,options,varargin);

Error in trial11 (line 70)

sol = bvp4c(@odefun, @bcfun, sol1,options);

Torsten 2024-2-23

编辑：Torsten 2024-2-23

在 MATLAB Online 中打开

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