Runge-Kutta 4th Order

I have t solve this equation (t^2)*y'' - 2*t*y' + 2*y = (t^3)*log(t) to solve first and secondly to compare the solutions with the theoretical solution y(t)= (7/4)*t + (t^3/2)*log(t) - (3/4)*t^3 ( 1<=t<=2, y'(1)=0 , y(1)=1, step h=0.05 ). I have a possible formula in my book but it doesn't seem to work. Is it possible for someone tο write an indicative formula, please?

[James Tursa]

What have you done so far? What specific problems are you having with your code? The phrase "... doesn't seem to work ..." doesn't really tell us much. Can you post your code?

[bk97]

>> syms y;
syms t;
format long;
f=inline('2*y'+((2*y)/t)+t^2*log(t)));
y1=1;
N=10;
y(1)=1;
for = i = 1:N
k1=h*f(y(i));
k2=h*f(y(i)+(k1)/2);
k3=h*f(y(i)+k3);
y(i+1)=y(i) + (k1 + 2*k2 + 2*k3 + k4)/6;
end

i adjusted my problem to the example of my book, please just help me with a new formula!

[James Tursa]

Why are you using symbolic variables? RK4 is a numerical technique used to solve initial value problems numerically. I.e., you would use it to find the values of your variables at future times. So I don't understand your use of symbols. What are you really trying to do here? E.g., seems like you could use the 2nd order ODE example in this link as a guide for solving your problem:

https://www.mathworks.com/help/matlab/ref/ode45.html?searchHighlight=ode45&s_tid=doc_srchtitle

[bk97]

 h=0.05;
 t=1:h:2;
 y=zeros(1,length(t));
 dy=diff(y,t);
 y(1)=1;
 dy(1)==0;
 f=@(t) 2*dy+(2*y)/t+t^2*log(t);
 for i=1:(length(t)-1)                             
    k1 = f(t(i),y(i));
    k2 = f(t(i)+0.5*h,y(i)+0.5*h*k1);
    k3 = f((t(i)+0.5*h),(y(i)+0.5*h*k2));
    k4 = f((t(i)+h),(y(i)+k3*h));

    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h; 
end

i have done this and this error apears "Error using diff Difference order N must be a positive integer scalar." can you help me?

Answer 1

[James Tursa]

OK, I will offer a bit more help here (well, actually a lot more help). Your most immediate problem is that you are treating your 2nd order ODE problem as if it is a 1st order ODE problem. So all of your stuff involving y(i) and y(i+1) etc is wrong because that is what you would do for a 1st order ODE (the result at each time step is a scalar). You need to look again at the 2nd order ODE example in this link that I have already given to you:

https://www.mathworks.com/help/matlab/ref/ode45.html?searchHighlight=ode45&s_tid=doc_srchtitle

In there you will see that for the 2nd order ODE problem the "y" solution at each step is in fact a 2-element vector. That was the hint that you needed to pull from this example. Your solution state is a 2-element vector, not a scalar. So all of your equations need to be rewritten with this in mind. The setup is as follows:

y is a 2-element vector
y(1) is your original y
y(2) is y'

So, start with your original equation:

(t^2)*y'' - 2*t*y' + 2*y = (t^3)*log(t)

Solve this for y''

y'' =  (2*t*y' - 2*y + (t^3)*log(t)) / t^2

or simplified a bit if you want

y'' =  2*y'/t - 2*y/t^2 + t*log(t)

So notionally if you start with

y(1) = y
y(2) = y'

Then the derivative of this (call it dy) will be

dy(1) = y'
dy(2) = y''

Substituting, you get

dy(1) = y(2)
dy(2) = 2*y(2)/t - 2*y(1)/t^2 + t*log(t)

These expressions for the elements of dy are the equivalent of what your "f" function is in your posted code. So you just need to rewrite your "f" function to take a 2-element y vector and time t, and output a 2-element dy vector. Then use that in your RK4 scheme, keeping in mind that at each point the state is a 2-element vector.

So, changing your scalar y to 2-element y and correcting some of your typos gives you the following framework:

   h = 0.05;
   t = 1:h:2;
   y = zeros(2,length(t));  % <-- changed 1 to 2, each column y(:,i) is the state at time t(i)
   y(1,1) = 1;  % <-- The initial value of y  at time 1
   y(2,1) = 0;  % <-- The initial value of y' at time 1
   f = ____________;  % <-- Derivative takes time t and a 2-element y and returns a 2-element result
   for i=1:(length(t)-1)  % At each step in the loop below, changed y(i) to y(:,i) to accomodate 2-element results
      k1 = f( t(i)      , y(:,i)         );
      k2 = f( t(i)+0.5*h, y(:,i)+0.5*h*k1);
      k3 = f( t(i)+0.5*h, y(:,i)+0.5*h*k2);
      k4 = f( t(i)+    h, y(:,i)+    h*k3);
      y(:,i+1) = y(:,i) + (1/6)*(k1 + 2*k2 + 2*k3 + k4)*h; 
   end 

Y = (7/4)*t + (t.^3/2).*log(t) - (3/4)*t.^3;  % <-- The analytic solution
figure
hold on
plot(t,Y,'b','LineWidth',5);  % Plot the analytic solution in thick blue
plot(t,y(1,:),'r','LineWidth',2);  % Plot the RK4 solution in thin red
grid on
legend('analytic','RK4');
xlabel('Time')
ylabel('y')
title('Solution to (t^2)*y'''' - 2*t*y'' + 2*y = (t^3)*log(t)');

The only thing you have left to do is fill in the code for the function f. I have given you what you need for this above.

[bk97]

you did a great job and i am very thankful for that.. all i need is to fill the gap with this? : 2*y(2)/t - 2*y(1)/t^2 + t*log(t) ?

[James Tursa]

Remember that the state of your system is a 2-element column vector y. The derivative function f takes a time t and a 2-element column vector y as inputs, and returns a 2-element column vector output. So you need to fill out this:

f = @(t,y) ________________________

The stuff you put in the blank should be a 2-element column vector derivative ... a function of t and y.

[bk97]

Well, Mr James i tried everything. Im getting really desperate and I can't find what the function in the blank space could be. Please, help me! Tell me exactly what the correct function is. Thank you in advance!

[James Tursa]

dy as a column vector ... remember you are passing in y and returning the derivative of y, which is the dy I spelled out for you above.

[Peng Li]

Define f as a common function (in another .m file) instead of an anonymous function, which is much easier to implement, as follows. Then the result for the code of @James Tursa can be obtained.

function dy = f( t,y )
  dy = zeros(2,1);
  dy(1) = y(2);
  dy(2) = 2*y(2)./t - 2*y(1)./(t.^2) + t.*log(t);
end

%

Answer 2

[Jan]

A Runge-Kutta as solver for a symbolic function? This is strange. Then your function depends on the inputs y and t, but inside your runge-Kutta-code you call it as f(x) only.

Start with transforming the 2nd order ODE to a set of equations in 1st order. Then omit the "syms", but create the solution numerically. You wil find many working examples when you search for "Matlab runge kutta".

[bk97]

yes i know this is wrong, but i'm working it for 6 hours and still can't do anything. If you have a code that works please show it to me and don't take into consideration my code! your help is important!

[Jan]

@bill: As long as I assume, that this is a homework, I do not like to post a working solution.

Do you know how to convert a higher order ODE to a set of 1st order ODEs? It is not complicated, or so to say trivial. Start with this. Then run a loop over the time t instead of "for = i = 1:N" (with an unintented = in between). Call your function as "f(y,t)". And as said before: There are many working Runge-Kutta codes in the net.

[bk97]

no i don't know this!

[bk97]

 h=0.05;
 t=1:h:2;
 y=zeros(1,length(t));
 dy=diff(y,t);
 y(1)=1;
 dy(1)==0;
 f=@(t) 2*dy+(2*y)/t+t^2*log(t);
 for i=1:(length(t)-1)                             
    k1 = f(t(i),y(i));
    k2 = f(t(i)+0.5*h,y(i)+0.5*h*k1);
    k3 = f((t(i)+0.5*h),(y(i)+0.5*h*k2));
    k4 = f((t(i)+h),(y(i)+k3*h));

    y(i+1) = y(i) + (1/6)*(k_1+2*k_2+2*k_3+k_4)*h; 
end

i have done this and this error apears "Error using diff Difference order N must be a positive integer scalar." can you help me?

[Jan]

@bill: The error message means, that in diff the 2nd input means the n'th difference, see doc diff. Therefore diff(y,t) is not meaningful here.

Start with defining a system of 1st order ODEs at first. See e.g. https://www.mathworks.com/help/symbolic/odetovectorfield.html or 2nd order ode to 1st order. Use a function instead of an anonymous function, because it might look easier.