How to solve runge kutta using implicit method

hey guys I am trying to write a script using runge kutta implicit method 4th order I got the following equations

for n =2:(N-1)  
k_1 = f(tout(n-1) +(0.5+(sqrt(3)/6)*h) , yout(:,n-1) + 0.25*h*k_1 + (0.25 + (sqrt(3)/6))*h * k_2);
k_2 = f(tout(n-1) +(0.5-(sqrt(3)/6)*h) , yout(:,n-1) + (0.25-(sqrt(3)/6))*h*k_1 + 0.25*h * k_2);
yout(:,n) = yout(:,n-1) + (h/2)*(k_1 + k_2); 
end

but to find k_1, I need a value for k_1 and to find k_2 I need a value for k_2

how I am supposed to do it? Didn't undestrand implict method...

Answer 1

[Sara]

Implicit means the equation has no analytic solution, i.e. you'll have to solve it iteratively. Meaning, you try guessing the value of your unknown, plug it into your equation and see if the right side is equal to the left side. In your case, for each iteration, you'll have to solve a system of algebraic equations. You can use fsolve for that, just rewrite your equations in k_1 and k_2 in the form f = (right side)-(left side)

[Rong]

Hi, thank you I tried to implement fsolve inside the loop and it didnt work do I need to create another script just to define these equations?

Answer 2

[Rong]

I did this

for n =2:(N-1)                              % calculation loop

      fun = @(k) [k(1) - (f(tout(n-1) +(0.5+(sqrt(3)/6)*h) , yout(:,n-1) + 0.25*h*k(1) + (0.25 + (sqrt(3)/6))*h * k(2)));
                        k(2) - f(tout(n-1) +(0.5-(sqrt(3)/6)*h) , yout(:,n-1) + (0.25-(sqrt(3)/6))*h*k(1) + 0.25*h * k(2))];  

      fsolve(fun, [h;h]);

      yout(:,n) = yout(:,n-1) + (h/2)*(k(1) + k(2)); % main equation    

end

and it says that

IRK4
Warning: Trust-region-dogleg algorithm of FSOLVE cannot handle non-square systems; using
Levenberg-Marquardt algorithm instead. 
> In fsolve at 314
  In IRK4Solver at 25
  In main at 38

No solution found.

fsolve stopped because the last step was ineffective. However, the vector of function
values is not near zero, as measured by the default value of the function tolerance.

[Sara]

What is the function f that you are trying to integrate? Can you attach all your code, input values included?