Почему числовая ошибка в моей реализации Runge-Kutta не уменьшается как N ^ a? - PullRequest
Новая
       12

Почему числовая ошибка в моей реализации Runge-Kutta не уменьшается как N ^ a?

2 голосов
/ 10 октября 2011

Я пытаюсь определить, сколько шагов требуется, чтобы метод Рунге-Кутта ("RK4") находился в пределах 0,01% от точного решения обыкновенного дифференциального уравнения. Я сравниваю это с методом Эйлера . Оба должны привести к прямой линии на графике журнала. Мое решение Эйлера кажется правильным, но я получаю кривое решение для РК. Они основаны на одном и том же коде, поэтому я полностью запутался в проблеме.

РЕДАКТИРОВАТЬ: Извините за удаление ссылок в Википедии. Это не позволило бы мне сохранить более одной ссылки, так как я новый пользователь. Однако оба метода подробно описаны в Википедии, как и моя реализация.

Если кто-то захочет заняться моей проблемой, код приведен ниже, а графики находятся в этом файле Word на dropbox.com . Да, это проблема с домашним заданием; Я публикую это, потому что я хотел бы понять, что не так в моем мыслительном процессе. 

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

Это мой код РК4: 

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);

Это мой код Эйлера: 

%# 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 голос
/ 10 октября 2011

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

Переменным, указанным в теле выражения [...], должно быть присвоено значение для них в то время вы создаете анонимную функцию, которая использует их. на конструкция, MATLAB фиксирует текущее значение для каждой переменной, указанной в тело этой функции. Функция будет продолжать связывать это значение с переменной, даже если значение должно измениться в рабочей области или выйти области видимости

Значение h в k1 ... k4 остается постоянным даже при изменении h в цикле while.

Одним из решений является добавление h к анонимным функциям: 

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