Second Order ODE with Runge Kutta 3 "K's problem"

Hi, I was given a matlab project on second order ODE.

Question: x''(1+sin(t))^2 = -cos(t), initail value x(0) = 2.

I changed it to first order and tried solving it using the following code , but the teacher said its wrong, that i need more k's, but i don't quite understand what the content of the k's are going to be. Please me check the code, Thanks in advance.

intmin=0;
intmax=pi/4;
h=0.2;
g=0.05;
numnodes=ceil(((intmax-intmin)/h)+1);
numnodes2=ceil(((intmax-intmin)/g)+1);

inival1 = 2;
inival2 = 0;

%0.2
t=zeros(1,numnodes);
x1=zeros(1,numnodes);
x2=zeros(1,numnodes);

t(1)=intmin;
x1(1)=inival1;
x2(1)=inival2;

f1 = @(t,x1,x2) x2;
f2 = @(t,x1,x2) (-cos(t))./((1+sin(t)).^2);

for i= 2:numnodes
    t(i) = t(i-1)+h;
    k1 = f2(t(i-1),x1(i-1),x2(i-1));
    %k11=;
    k2 = f2(t(i-1)+h/2,x1(i-1)+(h/2)*k1,x2(i-1)+(h/2)*k1);
    %k22=;
    k3 = f2(t(i-1)+h/2, x1(i-1)+(h/2)*k2,x2(i-1)+(h/2)*k2);
    %k33=;
    k4 = f2(t(i-1)+h, x1(i-1)+h*k3,x2(i-1)+h*k3);
    %k44=;
    x1(i) = x1(i-1)+(h/6)*(k1+2*k2+2*k3+k4);
    x2(i) = x2(i-1)+(h/6)*(k1+2*k2+2*k3+k4);
end

%0.05
o=zeros(1,numnodes2);
z1=zeros(1,numnodes2);
z2=zeros(1,numnodes2);

o(1)=intmin;
z1(1)=inival1;
z2(1)=inival2;

m1 = @(o,z1,z2) z2;
m2 = @(o,z1,z2) (-cos(o))./((1+sin(o)).^2);

for i= 2:numnodes2
    o(i) = o(i-1)+g;
    k1 = m2(o(i-1),z1(i-1),z2(i-1));
    %k11=;
    k2 = m2(o(i-1)+g/2,z1(i-1)+(g/2)*k1,z2(i-1)+(g/2)*k1);
    %k22=;
    k3 = m2(o(i-1)+g/2, z1(i-1)+(g/2)*k2,z2(i-1)+(g/2)*k2);
    %k33=;
    k4 = m2(o(i-1)+g, z1(i-1)+g*k3,z2(i-1)+g*k3);
    %k44=;
    z1(i) = z1(i-1)+(g/6)*(k1+2*k2+2*k3+k4);
    z2(i) = z2(i-1)+(g/6)*(k1+2*k2+2*k3+k4);
end

figure 
plot(t,x2,'.-',o,z2,'.-')
hold on

syms x(t)
dz = diff(x);

ode = diff(x,t,2) == (-cos(t))./((1+sin(t)).^2);
cond1 = x(0) == 2;
cond2 = dz(0) == 0;
dsolve(ode,cond1,cond2)

tt = linspace(0,1);
yy = 4 - 2./(tan(tt/2) + 1) - tt;

plot(tt,yy)
legend("Runge Kutta 3 0.2","0.05", "Actual");
hold off

Answer 1

[James Tursa]

You are using the same k's for x1 and x2, which is incorrect. That is, you are using f2( ) to calculate the k's, but these should only apply to the 2nd derivative, not to the 1st derivative. Your 1st derivative function f1( ) isn't even used, but it should be.

I would suggest you set up your state as a 2-element vector instead of separate variables x1 and x2. This will help you to write vectorized code for your RK4 scheme, and will also match what you would need to do when you move to the MATLAB function ode45( ). So, only have one derivative function, f( ), and have it work with a 2-element state vector. E.g.,

f = @(t,x) [ x(2); (-cos(t))./((1+sin(t)).^2) ];

Then use this to calculate your k's and to update your state at each step. Instead of individual scalars at each step, you will have 2-element vectors at each step. Maybe store them by columns. E.g.,

This

x1=zeros(1,numnodes);
x2=zeros(1,numnodes);

becomes this

x=zeros(2,numnodes); % 2-element state vectors stored by columns

And this

x1(1)=inival1;
x2(1)=inival2;

becomes this

x(:,1)=[inival1;inival2];

And this

k1 = f2(t(i-1),x1(i-1),x2(i-1));

becomes this

k1 = f(t(i-1),x(:,i-1));

And this

    x1(i) = x1(i-1)+(h/6)*(k1+2*k2+2*k3+k4);
    x2(i) = x2(i-1)+(h/6)*(k1+2*k2+2*k3+k4);

becomes this

    x(:,i) = x(:,i-1)+(h/6)*(k1+2*k2+2*k3+k4);

etc.

[Retr0]

Thanks James, I tried implementing your correction and I ended up having multiple errors, perhaps maybe if you can try copying the code and editing it yourself would be much helpful. Thank you very much again for the help :)

[James Tursa]

Post your updated code as a comment here and we can work on fixing the errors.

[Retr0]

intmin=0;
intmax=pi/4;
h=0.2;
g=0.05;
numnodes=ceil(((intmax-intmin)/h)+1);
numnodes2=ceil(((intmax-intmin)/g)+1);

inival1 = 2;
inival2 = 0;

%0.2
t=zeros(1,numnodes);
x=zeros(2,numnodes);

t(1)=intmin;
x(:,1)=[inival1;inival2];

f = @(t,x) [ x(2); (-cos(t))./((1+sin(t)).^2) ];

for i= 2:numnodes
    t(i) = t(i-1)+h;
    k1 = f(t(i-1),x(:,i-1));
    %k11=;
    k2 = f(t(i-1)+h/2,x(:,(i-1)+(h/2)*k1));
    %k22=;
    k3 = f(t(i-1)+h/2, x(:,(i-1)+(h/2)*k2));
    %k33=;
    k4 = f(t(i-1)+h, x(:,(i-1)+h*k3));
    %k44=;
    x(:,i) = x(:,i-1)+(h/6)*(k1+2*k2+2*k3+k4);
end

%0.05
o=zeros(1,numnodes2);
z=zeros(2,numnodes2);

o(1)=intmin;
z(:,1)=[inival1,inival2];

m = @(o,z) [ z(2); (-cos(o))./((1+sin(o)).^2)];

for i= 2:numnodes2
    o(i) = o(i-1)+g;
    k1 = m(o(i-1),z(:,i-1));
    %k11=;
    k2 = m(o(i-1)+g/2,z(:,(i-1)+(g/2)*k1));
    %k22=;
    k3 = m(o(i-1)+g/2, z(:,(i-1)+(g/2)*k2));
    %k33=;
    k4 = m(o(i-1)+g, z(:,(i-1)+g*k3));
    %k44=;
    z(:,i) = z(:,i-1)+(g/6)*(k1+2*k2+2*k3+k4);
end

figure 
plot(t,x2,'.-',o,z2,'.-')
hold on

syms x(t)
dz = diff(x);

ode = diff(x,t,2) == (-cos(t))./((1+sin(t)).^2);
cond1 = x(0) == 2;
cond2 = dz(0) == 0;
dsolve(ode,cond1,cond2)

tt = linspace(0,1);
yy = 4 - 2./(tan(tt/2) + 1) - tt;

plot(tt,yy)
legend("Runge Kutta 3 0.2","0.05", "Actual");
hold off

[James Tursa]

You don't have the indexing syntax correct. E.g., this line

k2 = f(t(i-1)+h/2,x(:,(i-1)+(h/2)*k1));

should be this

k2 = f( t(i-1) + h/2 , x(:,i-1) + (h/2)*k1 );

Spaces are your friends ... add them throughout your code to help readability.

[Retr0]

Its working now, but it gives the same graph as the first code. My teacher was said its wrong because we need to add more " K's ", that we should have 8 "K's" in the loop and I don't really understand what he meant by that.

[James Tursa]

You do have eight k's ... each of k1, k2, k3, k4 is a 2-element vector, so that makes eight k's total (i.e. four elements are for the 1st derivative and the other four elements are for the 2nd derivative). Please post your current updated code and let's work from that.

[Retr0]

ohh, I think I understand it now. Thanks, Here is the code.

intmin = 0;
intmax = pi/4;
h = 0.2;
g = 0.05;
numnodes = ceil(((intmax-intmin)/h) + 1);
numnodes2 = ceil(((intmax-intmin)/g) + 1);

inival1 = 2;
inival2 = 0;

%0.2
t = zeros(1, numnodes);
x = zeros(2, numnodes);

t(1) = intmin;
x(:,1) = [inival1; inival2];

f = @(t,x) [ x(2); (-cos(t))./((1+sin(t)).^2) ];

for i= 2:numnodes
    t(i) = t(i-1) + h;
    k1 = f(t(i-1), x(:,i-1));
    %k11=;
    k2 = f(t(i-1)+h/2, x(:,i-1) + (h/2)*k1);
    %k22=;
    k3 = f(t(i-1)+h/2, x(:,i-1) + (h/2)*k2);
    %k33=;
    k4 = f(t(i-1)+h, x(:,i-1) + h*k3);
    %k44=;
    x(:,i) = x(:,i-1) + (h/6)*(k1+2*k2+2*k3+k4);
end

%0.05
o = zeros(1,numnodes2);
z = zeros(2,numnodes2);

o(1) = intmin;
z(:,1) = [inival1,inival2];

m = @(o,z) [ z(2); (-cos(o))./((1+sin(o)).^2)];

for i= 2:numnodes2
    o(i) = o(i-1)+g;
    k1 = m(o(i-1), z(:,i-1));
    %k11=;
    k2 = m(o(i-1)+g/2, z(:,i-1) + (g/2)*k1);
    %k22=;
    k3 = m(o(i-1)+g/2, z(:,i-1) + (g/2)*k2);
    %k33=;
    k4 = m(o(i-1)+g, z(:,i-1) + g*k3);
    %k44=;
    z(:,i) = z(:,i-1) + (g/6)*(k1+2*k2+2*k3+k4);
end

figure 
plot(t,x,'.-',o,z,'.-')
hold on

syms x(t)
dz = diff(x);

ode = diff(x,t,2) == (-cos(t))./((1+sin(t)).^2);
cond1 = x(0) == 2;
cond2 = dz(0) == 0;
dsolve(ode,cond1,cond2)

tt = linspace(0,1);
yy = 4 - 2./(tan(tt/2) + 1) - tt;

plot(tt,yy)
legend("Runge Kutta 3 0.2","0.05", "Actual");
hold off

[James Tursa]

I need to leave and won't get back to this until later this evening. However, here is a slightly modified version of your code you can play with that has the MATLAB RK45 integrator function ode45( ) used and plotted next to your hand coded RK4 solution. Looks like a pretty good comparison to me. I would re-examine the DE you are using to make sure it is the same one you are supposed to be solving. And then also double check your symbolic solution to make sure it is also solving the same problem. Take the symbolic solution 2nd derivative and plug it into your DE to see that it actually solves the problem you think it is solving.

intmin = 0;
intmax = pi/4;
h = 0.2;
g = 0.05;
numnodes = ceil(((intmax-intmin)/h) + 1);
numnodes2 = ceil(((intmax-intmin)/g) + 1);
inival1 = 2;
inival2 = 0;
%0.2
t = zeros(1, numnodes);
x = zeros(2, numnodes);
t(1) = intmin;
x(:,1) = [inival1; inival2];
f = @(t,x) [ x(2); (-cos(t))./((1+sin(t)).^2) ];
for i= 2:numnodes
    t(i) = t(i-1) + h;
    k1 = f(t(i-1), x(:,i-1));
    %k11=;
    k2 = f(t(i-1)+h/2, x(:,i-1) + (h/2)*k1);
    %k22=;
    k3 = f(t(i-1)+h/2, x(:,i-1) + (h/2)*k2);
    %k33=;
    k4 = f(t(i-1)+h, x(:,i-1) + h*k3);
    %k44=;
    x(:,i) = x(:,i-1) + (h/6)*(k1+2*k2+2*k3+k4);
end
%0.05
o = zeros(1,numnodes2);
z = zeros(2,numnodes2);
o(1) = intmin;
z(:,1) = [inival1,inival2];
m = @(o,z) [ z(2); (-cos(o))./((1+sin(o)).^2)];
for i= 2:numnodes2
    o(i) = o(i-1)+g;
    k1 = m(o(i-1), z(:,i-1));
    %k11=;
    k2 = m(o(i-1)+g/2, z(:,i-1) + (g/2)*k1);
    %k22=;
    k3 = m(o(i-1)+g/2, z(:,i-1) + (g/2)*k2);
    %k33=;
    k4 = m(o(i-1)+g, z(:,i-1) + g*k3);
    %k44=;
    z(:,i) = z(:,i-1) + (g/6)*(k1+2*k2+2*k3+k4);
end
[TT,XX] = ode45(f,[intmin intmax],x(:,1));
figure 
plot(t,x(2,:),'.-',o,z(2,:),'.-',TT,XX(:,2),'.-')
hold on
syms x(t)
dz = diff(x);
ode = diff(x,t,2) == (-cos(t))./((1+sin(t)).^2);
cond1 = x(0) == 2;
cond2 = dz(0) == 0;
dsolve(ode,cond1,cond2)
tt = linspace(0,1);
yy = 4 - 2./(tan(tt/2) + 1) - tt;
%plot(tt,yy)
%legend('Runge Kutta 3 0.2','0.05', 'Actual');
legend('Runge Kutta 3 0.2','0.05','ode45');
grid on
hold off

[Retr0]

Ok, Thank you very much for the help. :)

[Retr0]

here is the complete question incase you need it.

Equation: x'' ( 1+sin(t) ) ^2 = - cos(t)

Initail value: x(0) = 2;

Method: Runge Kutta III

step sizes: h = 0.2, h=0.05

Thanks :)