Implementation implicit Gauss method

General Usage
,
1

I would like to implement Runge Kutta Gauss method. This method is implicit and need a nonlinear solver to do so. For that I define a fixed point solver FixIter, then the function Gauss(f::Function,u₀,h,N)

using LinearAlgebra
function FixIter(fun,u,tol)
    while true
        uold = u
        u = fun(u)
        norm(uold - u) > tol || break
    end
    return u
end
function  Gauss(f::Function,u₀,h,N)
A = [0.25 0.25-√3/6 ; 0.25+√3/6 0.25]

b = [0.5, 0.5] ;
u = zeros(length(u₀),N+1)
u[:,1]  =   u₀;
        
for i in 1:size(u,2)-1; 
 g(K) = [f(u[:,i]+ h * (A[1,1] * K[1] + A[1,2] * K[2])),f(u[:,i]+ h * (A[2,1] * K[1] + A[2,2] * K[2]))]
    K = FixIter(g,f(u[:,i]),1e-12);
 u[:,i+1] = u[:,i] + h * b[1] * K[1]  + h * b[2] * K[2];
end
    return u; 
end

My example is a pendulum problem with the vector field X = (y, \sin(x))

f(u) = [u[2], -sin(u[1])];
h=0.05; N=100;
x0=[1.,0.]
u = Gauss(f,x0,h,N)

However I get the following error

Screenshot from 2024-09-01 23-27-47
Screenshot from 2024-09-01 23-27-471317×508 82.6 KB

I appreciate any help you can provide.

2

if I’m understanding correctly, I believe this is just RadauIIA3 from OrdinaryDiffEq.jl

3

No, this is 2-stages Gauss method https://en.wikipedia.org/wiki/Gauss%E2%80%93Legendre_method

4

You are adding a vector to a scalar here?

PS. Please don’t post screenshots of text. Post quoted text.

6

I solve it but still not working

7

This code is also working but still give me a wrong result

function RKStep(f,x0,h,A,b)
 x=[] ; push!(x,x0);
tol = 1e-12;  
for i in 1:N-1; 
g(K) = f(x0 + h * A * K);
K = FixIter(g,x0,tol)
        x0 = x0 .+ h * b'* K;  
push!(x,x0);
end    
    return stack(x, dims=1)
end