# rung kutta 4th order

4 views (last 30 days)
Mj on 11 Nov 2020
Commented: Alan Stevens on 12 Nov 2020
Hi,
I have a code for a Runge Kutta 2th order analysis of a 2nd order ODE but it has not give me the right answer.
% d2a/dt2=(16*a^(1/2))/35 - 18/(35*a) - (6*a^3)/35
%initial condition: da/dt(x=-∞)=0 و a(x=0)=0.705
Here is my code below:
clear all
clc
h=0.2;
t=-2:h:0;
Nt=numel(t);
A=zeros(2,Nt);
A(:,1)=[0.705,0];
f = @(t, A) [A(2);(16*A(1).^(1./2))./35 - 18./(35.*A(1)) - (6.*A(1).^3)./35];
for i=1:Nt-1
k1 = h.*f(t(i), A(:,i));
k2 = h.*f(t(i) + 0.5.*h, A(:,i) + 0.5.*k1);
A(:,i+1) = A(:,i) + (1.0./6.0).*(k1 + 2.*k2 );
end
plot(t, A, '-o');
xlabel('x/hu')
ylabel('h/hu')
grid on

Alan Stevens on 11 Nov 2020
Your integration loop isn't that of a fourth order RK routine. Presumably, -2 is standing in for -infinity here, so you correctly have da/dt(t = -2) = 0. However, you also have a(t=-2) = 0.705, but you want a(t=0) = 0.705. The attached code achieves this ( I simply tried a few different values for a(t=-2) until a(t=0) became 0.705.
% d2a/dt2=(16*a^(1/2))/35 - 18/(35*a) - (6*a^3)/35
%initial condition: da/dt(x=-∞)=0 و a(x=0)=0.705
h=0.1;
t= -2:h:0;
Nt=numel(t);
A=zeros(2,Nt);
a0 = 1.9551; dadt0 = 0;
f = @(t, A) [A(2);(16*A(1).^(1./2) - 18./A(1) - 6.*A(1).^3)./35];
for i=1:Nt-1
k1 = f(t(i), A(:,i));
k2 = f(t(i) + 0.5.*h, A(:,i) + 0.5*h.*k1);
k3 = f(t(i) + 0.5.*h, A(:,i) + 0.5*h.*k2);
k4 = f(t(i) + h, A(:,i) + h*k3);
A(:,i+1) = A(:,i) + (h/6).*(k1 + 2.*k2 +2*k3 + k4);
end
disp(A(:,end))
plot(t, A, '-o',-2,0,'r*',0,0.705,'r*');
xlabel('x/hu')
ylabel('h/hu')
legend('a','da/dt')
grid on
Alan Stevens on 12 Nov 2020
Ok, here is the coding:
% d2a/dt2=(16*a^(1/2))/35 - 18/(35*a) - (6*a^3)/35
%initial condition: da/dt(x=-?)=0 ? a(x=0)=0.705
f = @(t, A) [A(2);(16*A(1).^(1./2) - 18./A(1) - 6.*A(1).^3)./35];
h=0.1;
t= -2:h:0;
a0 = 2.25; % Initial guess at a(t=-2)
a0 = fzero(@(a) Zfn(a, t, f), a0); % fzero passes a0 to Zfn and adjusts
% it until Zfn returns 0.
disp(a0)
IC = [a0 0];
A = RK4(IC,t,f); % Calculate A using final value of a0
disp(A(:,end))
figure
plot(t, A, '-o',-2,0,'r*',0,0.705,'r*');
xlabel('x/hu')
ylabel('h/hu')
legend('a','da/dt')
grid on
function Z = Zfn(a, t, f) % This function returns the difference between
% a(x=0) and 0.705
IC = [a, 0];
A = RK4(IC, t, f);
Z = A(1,end) - 0.705;
end
function A = RK4(IC, t, f) %
Nt = numel(t);
A=zeros(2,Nt);
A(:,1)=IC;
h = 0.1;
for i=1:Nt-1
k1 = f(t(i), A(:,i));
k2 = f(t(i) + 0.5.*h, A(:,i) + 0.5*h.*k1);
k3 = f(t(i) + 0.5.*h, A(:,i) + 0.5*h.*k2);
k4 = f(t(i) + h, A(:,i) + h*k3);
A(:,i+1) = A(:,i) + (h/6).*(k1 + 2.*k2 +2*k3 + k4);
end
end
Interestingly, there seem to be two possible initial values for a(t=-2) that give a(t=0) = 0.705. To find the other one set the initial guess for a0 to a0 = 1.