Написание функции для метода Implicit Runge-Kutta (четвертый порядок) - PullRequest
Новая
       63

Написание функции для метода Implicit Runge-Kutta (четвертый порядок)

5 голосов
/ 23 апреля 2020

Я пытаюсь составить функцию, которая решит систему ODES, используя неявный метод Рунге-Кутты (IRK) порядка 4, но у меня возникают проблемы с правильным определением моего l oop. Здесь мы определяем IRK как

enter image description here

Любой совет будет принят с благодарностью! 

function [tout,yout] = IRK4Solver(f,t,y0) 
t = t(:); % ensure that t is a column vector
N = length(t); 
h = (t(end)-t(1))/(N-1); % calculate h by assuming t gridpoints are equidistant
d = length(y0); % defines the size of y0
y0 = y0(:); % ensures that y0 is a column vector
y = zeros(d,N);% allocate the output array y
y(:,1) = y0; % assign y0 to the first column of y

% The entries of the following tableau are provided in the lecture notes
c = [1/2-sqrt(3)/6;
   1/2+sqrt(3)/6];
A = [1/4, 1/4-sqrt(3)/6;
     1/4+sqrt(3)/6, 1/4];
b = [1/2;1/2];

%calculate the loop
for n=1:N                           
    xi_1 = y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi_1)+h.*A(1,2)f(t(n)+c(2).*h,xi_2);
    xi_2 = y(:,n)+h.*A(2,1).*f(t(n)+c(1).*h,xi_1)+h.*A(2,2)f(t(n)+c(2).*h,xi_2);

    y(:,n+1) = y(:,n)+h.*b(1).*f(t(n)+c(1).*h,xi_1)+h.*b(2)f(t(n)+c(2).*h,xi_2);
end

tout = t;
yout = y;

ДАЛЬНЕЙШИЕ ПОПЫТКИ

Я включил команду fsolve в мою команду для l oop. И все же программа не будет работать. 

for n=1:N                           
 eq=@(xi) [xi(1:3)-(y(:,n)+h.*A(1,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(1,2)f(t(n)+c(2).*h,xi(1:3)));
     xi(1:3)-(y(:,n)+h.*A(2,1).*f(t(n)+c(1).*h,xi(1:3))+h.*A(2,2)f(t(n)+c(2).*h,xi(1:3)))];
     xistar=fsolve(eq,[1 1 1;1 1 1]);
    y(:,n+1) = y(:,n)+h.*b(1).*f(t(n)+c(1).*h,xistar(1:3))+h.*b(2)f(t(n)+c(2).*h,xistar(1:3));
end

Ответы [ 2 ]

0 голосов
/ 27 апреля 2020

Я использовал двухэтапную схему Хаммера и Холлингсворта (1955) четвертого порядка, полученную в «Р. Мей и Дж. Ное,« Численное решение обыкновенных дифференциальных уравнений: начальные задачи », Норт-Холл . Math. Stud., Vol. 83, pp. 1–94, Jan. 1984 ".

Мой код: 

function [tout,yout] = IRK4Solver(f,t,y0) 

t = t(:); % ensure that t is a column vector
N = length(t); 
h = t(2)-t(1); % calculate h by assuming t gridpoints are equidistant
d = length(y0); % defines the size of y0
y0 = y0(:); % ensures that y0 is a column vector
y = zeros(d,N);% allocate the output array y
y(:,1) = y0; % assign y0 to the first column of y

k1 = 0;
k2 = 0;

%calculate the loop
for n=2:N    
    k1 = f(t(n-1)+(0.5+sqrt(3)/6)*h,y(:,n-1)+1/4*k1+(1/4+sqrt(3)/6)*h*k2);
    k2 = f(t(n-1)+(0.5-sqrt(3)/6)*h,y(:,n-1)+1/4*k2+(1/4-sqrt(3)/6)*h*k1);
    y(:,n) = y(:,n-1)+(h/2)*(k1+k2);
end

tout = t;
yout = y;

Надеюсь, это то, что вы ищете.

0 голосов
/ 26 апреля 2020

Вам нужно включить команду fsolve внутри вашего для l oop

...