matlab微分方程

Question

我有下面的微分方程，我解不出来。

​

​

关于这个方程，我们知道以下几点：

D(r)是一个三阶多项式

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

D(2)=2D(1)

u(1)=450

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

其中K和Te是常数。

我想用一个矩阵来近似这个问题，然后我设法解决了类似的方程：

u(1)和u'(2)具有相同的极限条件。对于这个方程，我用中心差分近似u‘和u’，并使用r=1到r=2之间的有限差分法，然后把结果放在matlab中的矩阵A中，把极限条件放在matlab中的向量Y中，然后运行u=A\Y，得到u值是如何变化的。这是我的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

正如TroyHaskin在评论中指出的那样，一个人可以将D确定为一个常数因子，而这个常数因子无论如何都会在D'/D中抵消。换句话说:我们可以假设D(1)=1 (一个方便的数字)，因为D可以乘以任何常量。现在很容易找到系数(done with Wolfram Alpha)，多项式结果是

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

导数D'(r) = -6r^2+18r-12。(还有一种更智能的方法来查找多项式，从D‘开始，D’与已知的根是二次的。)

我可能会立即使用这些信息，计算一阶导数的系数k：

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碰巧已经存在，并且长度大于length(r)，则结果数组将具有无关的条目。(这是难以跟踪的错误的潜在来源。)