在Matlab中使用Runge-Kutta二阶打靶法的问题

Question

让我的RK2算法对某个二阶线性微分方程起作用，我遇到了一些问题。我已经在下面发布了我当前的代码(带有提供的参数)。由于某些原因，每次迭代时，y1的值与真实值的偏差会更大。我们将非常感谢您的任何意见。谢谢!代码：

    f = @(x,y1,y2)  [y2; (1+y2)/x];
    a = 1;
    b = 2;
    alpha = 0;
    beta = 1;
   n = 21;
   h = (b-a)/(n-1);
   yexact = @(x) 2*log(x)/log(2) - x +1;
   ye = yexact((a:h:b)');
   s = (beta - alpha)/(b - a);
   y0 = [alpha;s];
   [y1, y2] = RungeKuttaTwo2D(f, a, b, h, y0);
   error = abs(ye - y1);

   function [y1, y2] = RungeKuttaTwo2D(f, a, b, h, y0)
           n = floor((b-a)/h);  

           y1 = zeros(n+1,1); y2 = y1;
           y1(1) = y0(1); y2(1) = y0(2);

           for i=1:n-1
                   ti = a+(i-1)*h;
                   fvalue1 = f(ti,y1(i),y2(i));
                   k1 = h*fvalue1;

                   fvalue2 = f(ti+h/2,y1(i)+k1(1)/2,y2(i)+k1(2)/2);
                   k2 = h*fvalue2;

                   y1(i+1) = y1(i) + k2(1);
                   y2(i+1) = y2(i) + k2(2);
           end
    end

Answer 1

你的确切解决方案是错误的。你的微分方程有可能遗漏了一个减号。

y2'=(1+y2)/x的解决方案是y2(x)=C*x-1，然后是y1'=y2和y1(x)=0.5*C*x^2-x+D。

如果y2方程中的符号被反转，y2'=-(1+y2)/x，one将得到具有积分y1(x)=C*log(x)-x+D的y2(x)=C/x-1，其中包含给定的精确解。

0=y1(1) = -1+D  ==>  D=1
1=y1(2) = C*log(2)-1  ==  C=1/log(2)

此外，积分循环中的数组的长度为n+1，因此循环必须是从i=1到n。否则，最后一个元素仍然为零，这给出了第二个边界条件的错误残差。

纠正这一点并将计算扩大到一个割线步骤找到离散化的正确解决方案，因为常微分方程是线性的。对于步长为0.05的二阶方法，精确解的误差是由0.000285有界的。