D:/demona/2.0/manual/doc/2-31set/matlab/picard.m

Go to the documentation of this file.

clear
clc
n=65;
N=n*n;

%later with streamfunction

s=zeros(N,1);
v=zeros(N,1);
sold=zeros(N,1);
vold=zeros(N,1);

%fs(s,v,n)

%fv(s,v,n)

damping1=0.001;
damping2=1.0;
%jac=numJacV(s,v,n);


maxPicard=500;
maxNewton=10;

tol=1e-5;

%x=ones(N,1);
%ms(x,s,v,n)



for i=1:maxPicard
    
    disp('######################')
    disp('picard iter')
    i
    
    if mod(i,100)==0
        damping=damping1*1.1;
    end
    
    %start with vorticity
    
    for j=1:maxNewton
        
        f0=fv(sold,v,n);
        
        normF=norm(f0);
        
        
        if normF<tol
            disp('-----')
            disp('newton exit at iter')
            j-1
            break
        end
        
        deltax=zeros(N,1);
        
        %        deltax=gmres(@(x)mv(x,s,v,n),-f0,min(100,N),tol,min(1000,N));
        
        jac=numJacV(sold,v,n);
        
        deltax=jac\-f0;
        
        v=v+deltax;
        
        
    end
    v=damping1*v+(1-damping1)*vold;
    %plotP(v,n)
    %pause
    
    for j=1:maxNewton
        
        f0=fs(s,vold,n);
        
        normF=norm(f0);
        
        
        if normF<tol
            disp('-----')
            disp('newton exit at iter')
            j-1
            break
        end
        
        deltax=zeros(N,1);
        
        %        deltax=gmres(@(x)ms(x,s,v,n),-f0,min(100,N),tol,min(1000,N));
        
        jac=numJacS(s,vold,n);
        
        deltax=jac\-f0;
        
        s=s+deltax;
        
        
    end
    plotP(-s,n)
    max(max(abs(s)))
    pause(0.5)
    
    s=damping2*s+(1-damping2)*sold;
    
    
    
    if (norm(sold-s)<=tol)&&(norm(vold-v)<=tol)&&(i>1)
                disp('picard end')
        
        break
    end
    
    
    sold=s;
    vold=v;
   
    
    %gmres(@(x)ms(x,s,v,n),b,100,tol,maxit);
    
end

%v

---