Как я могу решить BVP с несколькими граничными условиями? - PullRequest
Новая
       54

Как я могу решить BVP с несколькими граничными условиями?

0 голосов
/ 03 апреля 2020

Я сталкиваюсь с некоторыми проблемами при решении моих дифференциальных уравнений методом прямой разности. У меня есть две области (от 0 до t и от t до 1, а t это любое значение от 0 до 1). В первой области (от 0 до t) у меня есть 4 дифференциальных уравнения, а во второй области (от t до 1) у меня есть 2 дифференциальных уравнения, всего у меня шесть граничных условий (B C), 2 B C в первая область, 4 B C на границе раздела (t) и 1 B C во второй области.

Я следовал документации в MathWorks, чтобы решить BVP с несколькими граничными условиями (https://www.mathworks.com/help/matlab/math/solve-bvp-with-multiple-boundary-conditions.html)

Когда я запускаю свой код, я получаю сообщение об ошибке, недостаточно входные аргументы. Я думаю, что проблема заключается в регионе. Решатель BVP не получает необходимую область для решения уравнений. Пожалуйста, кто-нибудь, подскажите, как преодолеть эту проблему? это будет оценено. 

t=0.9;    
x = [0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 t t 0.95 1];
solinit = bvpinit(x,[zeros(1,90)]);
sol = bvp5c(@ex1ode,@ex1bc,solinit);
y = deval(sol,x);
plot(sol.x,sol.y(1,:));
function dydx = ex1ode(x,y,region)
h = 0.0001; t = 0.9; k = 100; Bi = 0.01; Da = 0.0001;
A0 = 100; A1 = 0.944545454545455; A2 = -0.4445454545454545; A3 = 7.3001203377370034e-43; A6 = 0.0101;
E1 = 0.002414161912718531; E2 = 3.55998187608043e-45 ; E3 = 0.00476100894331544; E4 = 0; E5 = -0.816974047611677;
D1 = 6.785528299098551; D2 = -10.83936751884141; D3 = 7.676961726023331; D4 = -2.03191960987143; D5 = -16.67685194071544;
uc=((-(x^2)/2)+A1*x+A2)/(-0.166667+(0.5*A1)+A2-(A2*t)+(Da*t)-(0.5*A1*(t^2))+0.166667*(t^3)+ 
(A3*sinh(A0*t)/A0));
up =(A0*(A3*cosh(A0*x)+Da))/(A0(-0.166667+A2-(A2*t)+(Da*t)+0.166667*(t^3)+A1*(0.5-0.5* 
(t^2)))+A3*sinh(A0*t));
yini=-E1*(x^2)-E2*cosh(A0*x)-E3*cosh(sqrt(A6)*x)-E4*cosh(2*A0*x)-E5;
tini=-D1*x-D2*(x^2)-D3*(x^3)-D4*(x^4)-D5;
n=90;
dydx=zeros(6,1);
for i=1:6:n
    switch region
        case 1  % x in [0 t]
            dydx(i)=y(i+1);
            dydx(i+1)=up*((y(i)-yini)/h)-(Bi./k)*(y(i+2)-y(i));
            dydx(i+2)=y(i+3);
            dydx(i+3)=Bi*(y(i+2)-y(i));
            yini=y(i);
        case 2  % x in [t 1]
            dydx(i+4)=y(i+5);
            dydx(i+5)=uc*((y(i+4)-tini)/h);
            tini=y(i+4);
end
end
end

function res=ex1bc(YL,YR)
n=90;k =0.01;
res=zeros(n,1);
for i=1:6:n
    res(i)=YL(i+1,1);
    res(i+1)=YL(i+3,1);
    res(i+2)=YR(i,1)-YR(i+2,1);
    res(i+3)=YR(i+2,1)-YL(i+4,2);
    res(i+4)=k*YL(i+5,2)-k*YR(i+1,1)-YR(i+3,1);
    res(i+5)=YR(i+5,2)-1;
end
end
...