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

### Retr0 (view profile)

on 15 Apr 2019
Latest activity Commented on by Retr0

on 16 Apr 2019

### James Tursa (view profile)

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

### James Tursa (view profile)

on 15 Apr 2019

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.

James Tursa

### James Tursa (view profile)

on 15 Apr 2019
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

### Retr0 (view profile)

on 15 Apr 2019
Ok, Thank you very much for the help. :)
Retr0

### Retr0 (view profile)

on 16 Apr 2019
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 :)