WIP Sg fg complete2 d

PDE/fokkerplanck2_complete.m

@@ -42,7 +42,7 @@
 
 %%
 % Select 6.1, 6.2, 6.3, etc where it goes as 6.test
-test = '6p1b'; 
+test = 'full'; 
 
 %%
 % Define a few relevant functions
@@ -56,7 +56,7 @@
         E = 0.0025;
         tau = 10^5;
     case '6p1b'
-        delta = 0.042;
+        delta = 3;
         Z = 1;
         E = 0.25;
         tau = 10^5;
@@ -65,6 +65,11 @@
         Z = 1;
         E = 0.0025;
         tau = 10^5;
+    case 'full'
+        delta = 0.3;
+        Z = 5;
+        E = 0.4;
+        tau = 10^5;
 end
 gamma = @(p)sqrt(1+(delta*p).^2);
 vx = @(p)1/vT*(p./gamma(p));
@@ -107,6 +112,8 @@
                     ret = ret + h(l) * P;
 
                 end
+            case 'full'
+                ret = x.*0 + 1;
         end
     end
 
@@ -129,12 +136,23 @@
                 ret = 2/(sqrt(pi)*a^3) * exp(-x.^2/a^2);
             case '6p1c'
                 ret = 2/(3*sqrt(pi)) * exp(-x.^2);
-
+            case 'full'
+                N = 1000;
+                h = 20/N;
+                Q = 0;
+                Fun = @(p)exp(-2/delta^2*sqrt(1+delta^2*p.^2));
+                for i = 1:N
+                    x0 = (i-1)*h;
+                    x1 = i*h;
+                    [xi,w] = lgwt(20,x0,x1);
+                    Q = Q+sum(w.*Fun(xi).*xi.^2);
+                 end
+                ret = exp(-2/delta^2*sqrt(1+delta^2*x.^2))/(2*Q);
         end
     end
 
 
-    function ret = soln_z(x,t)
+    function ret = soln_z(x,~)
         ret = x.*0+1;
     end
     function ret = soln_p(x,t)
@@ -143,8 +161,8 @@
 
 %% Setup the dimensions 
 
-dim_p.domainMin = 0.1;
-dim_p.domainMax = +10;
+dim_p.domainMin = 0;
+dim_p.domainMax = +20;
 dim_p.init_cond_fn = @(x,p,t) f0_p(x);
 
 dim_z.domainMin = -1;
@@ -236,7 +254,7 @@
 termE1_z = TERM_1D({pterm1});
 
 pterm1 = MASS(g2); 
-pterm2 = GRAD(num_dims,g3,1,'N','D');% Lin's Setting
+pterm2 = GRAD(num_dims,g3,1,'N','N');% Lin's Setting
 termE1_p = TERM_1D({pterm1,pterm2});
 
 termE1 = TERM_ND(num_dims,{termE1_p,termE1_z});
@@ -270,7 +288,7 @@
 g3 = @(x,p,t,dat) 1-x.^2;
 
 pterm1   = MASS(g1);% This is not needed - by Lin
-pterm2   = GRAD(num_dims,g2,1,'N','D');% Lin's Setting
+pterm2   = GRAD(num_dims,g2,1,'N','N');% Lin's Setting
 
 termR1_p = TERM_1D({pterm1,pterm2});
 
@@ -300,7 +318,7 @@
 %%
 % Add terms to the pde object
 
-pde.terms = {termC1,termC2,termC3,termE1,termE2,termR1,termR2};
+pde.terms = {termC1, termC2, termC3, termE1, termE2, termR1, termR2};
 
 %% Construct some parameters and add to pde object.
 %  These might be used within the various functions below.