Решение кусочно заданой системы дифференциальных уравнений - Matlab Опции

[X-Shadow]

Здравствуйте. Столкнулся с проблемой при решении системы дифференциальных уравнений которая представлена на рисунке в приложении. Эта система содержит два интервала. Из за того что это система уравнений второго порядка и я не знаю значений для первых производных на границах (y'(-20) и y'(10)) то это задача с граничными условиями и в теории может быть решена с помощью функции bvp4c. Проблема в том как прописать эту систему уравнений вместе с условием непрерывности в районе точки 0 (снизу изображения). Я попытался сделать следующее (функция для передачи в bvp4c):

Код

function [ dydx ] = tode( x , y )
dydx = zeros(2,1); % create zero array
dydx(1) = y(2);
if (x > 0)
    dydx(2) = b/a*((1+y(1))^(3/2)-1);
elseif (x <= 0)
    dydx(2) = 0;
end;
end

Однако эта функция не включает в себя условие с первыми производными. Подскажите пожалуйста как можно решить данную проблему? Спасибо.

Эскизы прикрепленных изображений

 

[X-Shadow]

Большое спасибо за помощь!

Потенциал на границе х=10 должен быть равен нулю т.к там металл который подключен к земле. Производная тоже должна быть равна нулю т.к внутри металла электрического поля нет (рассматриваем идеальный случай).

Касательно почему не в пакетах TCAD. Я сначала хочу получить аналитический результат а потом сравнить его с симуляцией.

Сейчас пытаюсь реализовать это в MATLAB но пока как то туго. Собственно код:

Код

global coef_A coef_B;

coef_A = 1.5477e+10;
coef_B = 9.7420e+03;

U = 1;              %значение y(-20e-7)
dy1 = 0;           %значение y'(-20e7)

err = 1e-1;       %допустимая ошибка значения y(10e-7)
end_val = 2;     %значение в точке y(10e-7)

xspan_neg = [-20e-7 : 1e-8 : -1e-8];      %интервалы
xspan_pos = [1e-8 : 1e-8 : 10e-7];

while end_val > err         %пока конечное значение больше ошибки
    
    y0_neg = [U dy1];                                                         %вектор потенциала и производной в точке -20у-7
    [X1,Y1] = ode45('neg_diff_eq',xspan_neg,y0_neg);          %вычисляем функцию y(x) и y'(x) на участке [-20e-7; 0]

    U1 = Y1(end,1);                                                             %значение y(+0), первое условие непрерывности
    dU1 = 0.5391*Y1(end,2);                                                %значение y'(+0), второй условие непрерывности

    y0_pos = [U1 dU1];                                                        %формируем вектор
    [X2,Y2] = ode45('pos_diff_eq',xspan_pos,y0_pos);           %вычисляем функцию y(x) и y'(x) на участке [0; 10e-7]

    end_val = abs(Y2(end,1));                                              %полученное значение значение y(10e-7)
    dy1 = dy1 - 1e4;                                                            %уменьшаем значение первой производной
    
end


X = [X1; X2]; Y = [Y1; Y2];                                                %склеиваем графики

subplot(2,1,1);                                                                  %рисуем
plot(X,Y(:,1));
subplot(2,1,2);
plot(X,Y(:,2));

А сами функции выглядят так:

Код

function [ dydx ] = neg_diff_eq( x , y )

dydx = zeros(2,1); % create zero array
dydx(1) = y(2);
dydx(2) = 0;
            
end

function [ dydx ] = pos_diff_eq( x , y )

global coef_A coef_B;

dydx = zeros(2,1); % create zero array
dydx(1) = y(2);

dydx(2) = coef_A*(1+coef_B*y(1))^(3/2)-1;
            
end

Но в итоге решение почемуто не сходится. Матлаб просто зависает пытаясь это решить

Сообщение отредактировал X-Shadow - Jan 17 2014, 03:23