Asked by Retr0
on 15 Apr 2019

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

Accepted Answer

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

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

Retr0
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 :)

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 0 Comments

Sign in to comment.