1
$\begingroup$

I have written the numerical scheme of PSOR (Projected Successive Over Relaxation) method of a obstacle program in 1D, which is as follows:

\begin{align*} &u_0=u_N=0\\ &u_i^{k+\frac{1}{2}}=(1-\omega)u_i^k+\omega\Big[ \frac{1}{2}\Big(u_{i-1}^{k+1}+u_{i+1}^k\big)+\frac{h^2}{2T}f_i\Big]\\ &u_i^{k+1}=\min\{g_i,u_i^{k+\frac{1}{2}}\} \end{align*}

where $g_i$ represents the obstacle function in $x_i$. For this case we can take $g$ as $$g=0.1(x-1)^2+0.02(x-1)+\frac{0.1}{(1+30(x-1)^2)}$$

I have tried to code it in Mathematica, however the output it's not correct. Could someone try to identify what is wrong with it? Thank you in advance!

    MetPSOR1D[a_,b_,n_,Omega_,t_,tol_,itMax_] := (h=(b-a)/n;
    
    T = Table[a+i*h,{i,0,n}];
        
    With[{list = Table[0,{n+1}]},
        Compile[{}, 
            Module[{listold = list, listnew = list, error = tol + 1, f = 1},
                
                For[k = 0, error > tol && k < itMax, k++, 
                    listnew[[1]]=0;
                    listnew[[n+1]]=0;
                    
                    For[i = 2, i < n+1, i++, 
                        g = 0.1*(T[[i-1]]-1)^2 + 0.02*(T[[i-1]]-1) + 0.1/(1 + 30*(T[[i-1]]-1)^2);
                        u12 = (1-Omega)listold[[i]] + (Omega(listold[[i-1]] + listold[[i+1]] + (f/t) * h^2))/2; 
                        aux = Min[{u12,g}];
                        listnew[[i]] = aux;
                    ]

                    error = Max[Abs[listnew - listold]];
                    listold = listnew;
                ]; 
                s = listnew;
            ]
        ]
    ]
    ); s
Improve this question
$\endgroup$
7
  • $\begingroup$ Could you show us what you’ve tried so far, in Mathematica code? $\endgroup$
    – MarcoB
    Dec 24, 2020 at 1:06
  • $\begingroup$ Yes, just uploaded it. Thank you for your help in advance! $\endgroup$
    – Alice
    Dec 24, 2020 at 1:40
  • 4
    $\begingroup$ There seems to be a critical semicolon missing between the assignment to T and the beginning of your Module expression. $\endgroup$
    – m_goldberg
    Dec 24, 2020 at 2:35
  • 3
    $\begingroup$ Well, you will have update your question with the corrected code and explain what your main issue is. Otherwise, this question will most likely be closed without an answer. $\endgroup$
    – m_goldberg
    Dec 24, 2020 at 10:58
  • 1
    $\begingroup$ Your code still seems to be missing a critical semicolon between the For[...] and error = .... Also, you need to supply us with test data so we can run your code. $\endgroup$
    – m_goldberg
    Dec 29, 2020 at 15:09

0

Reset to default