rung kutta 4th order

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

Answer 1

[Alan Stevens]

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;
A(:,1)=[a0, dadt0];

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

[Mj]

thanks alot. but i dont have the value of a(t=-2) and it is not equal to 0.705.

i just have the value of a(t=0)=0.705 and da/dt(t = -2) = 0

[Alan Stevens]

Exactly! But with the value of a(t=-2) set to 1.9551 you end up with the value a(t=0) of 0.705.

[Mj]

how do you underestand to set the value of a(t=-2) to 1.9551?

[Mj]

isn it possible to solve the equation without the value of a(t=-2)???

[Alan Stevens]

As I said, I simply tried a few values for a(t=-2), It was fairly easy to home in on the value of 1.9551.

To do this more formally you would select a single guess value for a(t=-2) and then use a search algorithm, such as Matlab's fzero function, to home in on an accurate value. However, it was so easy to find 1.9551 in this case that I didn't think it was worth writing the extra coding to make use of fzero. However, I recommend that you look up the documentation on fzero.

[Mj]

thanks.is it possible for you to write the code of using fzero for getting the value of 1.9551?

i am beginner in matlab.

[Alan Stevens]

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.