дифференциальное уравнение matlab

Question

У меня есть следующее дифференциальное уравнение, которое я не могу решить.

Мы знаем следующее о уравнением:

D (г) представляет собой третий класс Полином

D '(1) = D' (2) = 0

D (2) = 2D (1)

U (1) = 450

и '(2) = - K * (u (2) -Te)

Где K и Te - постоянные.

Я хочу приблизить проблему, используя матрицу, и мне удалось решить аналогичное уравнение: с теми же предельными условиями для u (1) и u '(2). В этом уравнении я аппроксимировал u 'и u' 'с центральными отличиями и использовал метод конечных разностей между r = 1 и r = 2. Затем я поместил результаты в матрицу A в matlab и предельные условия в векторе Y в matlab и запустил u = A \ Y, чтобы узнать, как изменяется значение u. Heres мой код MATLAB уравнения мне удалось решить:

clear 
    a=1; 
    b=2; 
    N=100; 
    h = (b-a)/N; 
    K=3.20; 
    Ti=450; 
    Te=20; 

    A = zeros(N+2); 
    A(1,1)=1; 
    A(end,end)=1/(2*h*K); 
    A(end,end-1)=1; 
    A(end,end-2)=-1/(2*h*K); 


    r=a+h:h:b; 

    %y(i) 
    for i=1:1:length(r) 
     yi(i)=-r(i)*(2/(h^2)); 
    end 
    A(2:end-1,2:end-1)=A(2:end-1,2:end-1)+diag(yi); 

    %y(i-1) 
    for i=1:1:length(r)-1 
     ymin(i)=r(i+1)*(1/(h^2))-1/(2*h); 
    end 
    A(3:end-1,2:end-2) = A(3:end-1,2:end-2)+diag(ymin); 

    %y(i+1) 
    for i=1:1:length(r) 
     ymax(i)=r(i)*(1/(h^2))+1/(2*h); 
    end 
    A(2:end-1,3:end)=A(2:end-1,3:end)+diag(ymax); 


    Y=zeros(N+2,1); 
    Y(1) =Ti; 
    Y(2)=-(Ti*(r(1)/(h^2)-(1/(2*h)))); 
    Y(end) = Te; 
    r=[1,r]; 
    u=A\Y; 
    plot(r,u(1:end-1)); 

Мой вопрос, как я могу решить дифференциальное уравнение первого?

Answer 1

Как указал Трой Хаскин в комментариях, можно определить D до постоянного множителя, и постоянный коэффициент в любом случае отменяется в D '/ D. Иными словами, можно считать, что D (1) = 1 (удобное число), так как D можно умножить на любую константу. Теперь легко найти коэффициенты (done with Wolfram Alpha), а полином оказывается

D (г) = 2г^3 + 9R^2-12r + 6

с производной D '(г) = -6r^2 + 18r-12. (Существует также более разумный способ найти многочлен, начиная с D», который является квадратичным известными корнями.)

я бы, вероятно, использовать эту информацию сразу, вычисляя коэффициент к первой производной:

r = a+h:h:b; 
k = 1+r.*(-6*r.^2+18*r-12)./(-2*r.^3+9*r.^2-12*r+6); 

Кажется, что k всегда положительно на интервале [1,2], поэтому, если вы хотите свести к минимуму изменения существующего кода, просто замените r (i) на r (i)/k (i) в нем ,

---

Кстати, вместо петель, как

for i=1:1:length(r) 
    yi(i)=-r(i)*(2/(h^2)); 
end 

один обычно делает просто

yi=-r*(2/(h^2)); 

Это векторизации делает код более компактным и может принести пользу производительность слишком (не так много в вашем примере, где решение линейной системы является узким местом). Другим преимуществом является то, что yi правильно инициализирован, а при построении цикла, если yi уже существует и имеет длину больше длины (r), результирующий массив будет иметь посторонние записи. (Это потенциальный источник труднодоступных ошибок.)