为什么龙格-库塔的数字错误在我实现不减少像N ^ ?

I'm trying to determine how many steps it takes for the Runge-Kutta Method ("RK4") to be within 0.01% of the exact solution of an ordinary differential equation. I'm comparing this to the Euler method. Both should result in a straight line on a loglog plot. My Euler solution seems to be correct but I am getting a curved solution for RK. They are based on the same code so I am completely confused about the problem.

我试图确定龙格-库塔法(“RK4”)在一个常微分方程的精确解的0.01%范围内需要多少步。我把它和欧拉法比较。两者都应该在loglog图上形成一条直线。我的欧拉解决方案似乎是正确的，但我正在为RK得到一个弯曲的解决方案。它们基于相同的代码，所以我完全不理解这个问题。

EDIT: Sorry for removing the wikipedia links. It wouldn't let me keep more than one link since I'm a new user. However, both methods are detailed on Wikipedia similarly to my implementation.

编辑:抱歉删除*链接。它不会让我保存多个链接，因为我是一个新用户。然而，这两种方法在*上都是类似于我的实现。

Should someone wish to tackle my problem, the code is below and graphs are in this Word file on dropbox.com. Yes this is a homework problem; I'm posting this because I'd like to understand what is wrong in my thought process.

如果有人想解决我的问题，下面的代码和图在这个Word文件的dropbox.com。是的，这是一个家庭作业问题;我发布这篇文章是因为我想知道我的思考过程中有什么错误。

f = @(x,y) x+y; %this is the eqn (the part after the @(t,y)

This is my RK4 code:

这是我的RK4代码:

k1=@(x,y) h*f(x,y);
k2=@(x,y) h*f(x+1/2*h,y+1/2*k1(x,y));
k3=@(x,y) h*f(x+1/2*h,y+1/2*k2(x,y));
k4=@(x,y) h*f(x+h,y+k3(x,y));

clear y x exact i
x(1)=0;
y(1)=2;
xn=1;
x0=0;

exact=3*exp(xn)-xn-1; %exact solution at x=1
%# Evaluate RK4 with 1 step for x=0...1
N=1; %# number of steps
h=(xn-x0)/N; %# step size
i=1;
y(i+1)=y(i)+1/6*k1(x(i),y(i))+1/3*k2(x(i),y(i))+1/3*k3(x(i),y(i))+1/6*k4(x(i),y(i));
RK4(N)=y(i+1);  %# store result of RK4 in vector RK4(# of steps)
E_RK4(N)=-(RK4(N)-exact)/exact*100;%keep track of %error for each N
Nsteps_RK4(N)=N;


%# repeat for increasing N until error is less than 0.01%
while -(RK4(N)-exact)/exact > 0.0001
    i=1;
    N=N+1;
    h=(xn-x0)/N;
    for i=1:N
        y(i+1)=y(i)+1/6*k1(x(i),y(i))+1/3*k2(x(i),y(i))+1/3*k3(x(i),y(i))+1/6*k4(x(i),y(i));
        x(i+1)=x(i)+h;
    end
    RK4(N)=y(i+1);
    Nsteps_RK4(N)=N;
    E_RK4(N)=-(RK4(N)-exact)/exact*100; %# keep track of %error for each N
end

Nsteps_RK4(end);

This is my Euler code:

这是我的欧拉代码:

%# Evaluate Euler with 1 step for x=0...1
clear y x i
x(1)=0;
y(1)=2;
N=1; %# number of steps
h=(xn-x0)/N; %# step size
i=1;
y(i+1)= y(i)+h*f(x(i),y(i));
Euler(N)=y(i+1); %# store result of Euler in vector Euler(# of steps)
E_Euler(N)=-(Euler(N)-exact)/exact*100;%# keep track of %error for each N
Nsteps_Euler(N)=N;
%# repeat for increasing N until error is less than 0.01%
while -(Euler(N)-exact)/exact > 0.0001
    i=1;
    N=N+1;
    h=(xn-x0)/N;
    for i=1:N
        y(i+1)= y(i)+h*f(x(i),y(i));   
        x(i+1)=x(i)+h;
    end
    Euler(N)=y(i+1);
    Nsteps_Euler(N)=N;
    E_Euler(N)=-(Euler(N)-exact)/exact*100; %# keep track of %error for each N
end

1 个解决方案

#1

See: http://www.mathworks.com/help/techdoc/matlab_prog/f4-70115.html#f4-71621

参见:http://www.mathworks.com/help/techdoc/matlab_prog/f4 # f4 - 71621 - 70115. - html

Variables specified in the body of the expression [...] must have a value assigned to them at the time you construct an anonymous function that uses them. Upon construction, MATLAB captures the current value for each variable specified in the body of that function. The function will continue to associate this value with the variable even if the value should change in the workspace or go out of scope.

在表达式的正文中指定的变量[…在构建使用它们的匿名函数时，必须赋予它们一个值。在构建过程中，MATLAB获取了该函数体中指定的每个变量的当前值。函数将继续将该值与变量关联，即使该值应该在工作区中更改或超出范围。

The value of h in k1...k4 remains constant even though you change h in the while loop.

h在k1中的值…k4保持不变，即使你在while循环中改变h。

One solution is to add h to the anonymous functions:

一个解决方案是将h添加到匿名函数中:

k1=@(x,y,h) h*f(x,y);

#1