Newton-galerkin method for non-linear differential equation

Hi

I want to implement the newton galerking method in order to solve the equation c(0)=1 & .I know there might be a built in function in matlab but I want to build my own code. The problem is that I am stuck in discretization of the equation for the R and J matrix(J δu=-R). I analyze c in quadratic basis functions , c=. The R must be nnx x 1 in dimensions in order to solve the system with newton iterative method but if I substitute c the R discretization has terms with i and also j so the matrix should be nnx x nnx in dimensions right?But thats impossible cause the R must have only i in discretization. I was used to implementing the galerkin method in linear differential equations where we solve the system Rij=0 and there is no need for initial values and iterative methods but now the R must be a vector not a matrix and I dont know how to deal with that. I know my question is a result of not fully undestanding the method. Any answer would help.

DISCRETIZATION: Ri=

J=

[J. Alex Lee]

Sounds like in the linear case you got some vocab and denotations a bit confused...in G/FEM you don't discretize equations; what you have shown as "DISCRETIZATION:..." should just be calculus if you have done it right.

You are discretizing your solution using basis functions. You've indexed them as "j" within an element, so actually your elemental jacobian should be dRj/dci, unless by convention you want them flipped (which could be the case). So it might just be a matter of clarifying how you want everything to be indexed and keeping the indices straight?

In the linear case, you were not letting R be a matrix, that wouldn't make sense. You had no need for R at all. It would just have been the coefficients of your linear model, which would be teh same as dRj/dci if you had indexed everything correctly.

For each finite element, you would have a list (vector) of all c_j (the solution vector), and an associated residual vector, with number of entries equal to the number of basis functions (and therefore number of elemental nodes) NNE, and Jacobian is a matrix of size NNExNNE. At each iteration you are just solving the "linearized" version of your actual problem at the current solution vector guess.

Does that help? It's been a while since I've gone through all that myself, so maybe I'm not being super clear...

[christos d]

Yeah ok i get it i think i disagree only in the part that it should be dRj/dci instead of dRi/dcj which is not the case since I discretized my solution with j, so the derivative should be with respect to j. Nevertheless, thank yor for your time spent replying to me and your help. I solved the problem.The issue on the code was that in dRi/dcj and Ri the initial guess included in them was calculated by mistake with local numbering so since i was using quadratic basis functions the program actually used only the first 3 values of the last approximation of cj no matter what the element was. Again thanks.

[J. Alex Lee]

Sure thing. It really doesn't matter what you call the index as long as you are using consistent indices throughout. It is less confusing to Express your discretized solution as a sum over $i$ if you want to index your solution and residual by $i$, but again you can call it whatever you want.

[christos d]

Hello again,

now that I found the solution of the differential equation in each node, do you have any idea how would I find the derivative of the solution using the basis functions.Obviously I can calculate it using finite differences easiily but thats not what I am looking for.I should utilize the finite element. I know that but how do implement it in code?

[J. Alex Lee]

It's based on the gradients of your basis functions, that's the whole idea behind the method. The algorithm anyway needs to compute gradients that define the ode/pde, so quickest way would be to just export those intermediate results if possible. If you need to do it outside of your solver I think it's as easy as evaluating the gradients of your basis functions and applying your solution weights as usual...but it's been so long I forgot details...

[christos d]

Yeahp it worked perfectly mr J Alex Lee. I calculated the derivative in each element somewhere in between(ξ=0.25 for input in basis functions) ,cause at the start or at the end of each element basis functions are not differentiable, using the solutions I got and the results are pretty much values that i expected so it must be right. It appears that derivative calculation is an easy task as soon as you create the code for the solution of equation.Thanks again for your help sir!

[J. Alex Lee]

My pleasure, glad it worked!