Comparison with Stefano's code for different meshs. Finer meshing for safety.

Code can be run for nu=10^-15 or 10^-6, significant deviations for precision=1.

@Ternary_model/Ternary_model.m

@@ -9,6 +9,7 @@ classdef Ternary_model < handle
         g0
         geometry = 2; % 2 ... Radial, 1 ... Cylindrical, 0 ... Slab
         ic = 'FRAP'
+        precision = 5;
         sol
         t
         tau
@@ -25,10 +26,11 @@ classdef Ternary_model < handle
         x0 = 1; % Center of Gauss initial condition
     end
     methods
-        function T = Ternary_model(geometry, ic, params, t)
+        function T = Ternary_model(geometry, ic, params, t, prec)
             T.geometry = geometry;
             T.ic = ic;
             T.t = t;
+            T.precision = prec;
             if ~isempty(params)
                 T.a = params(1);
                 T.b = params(2);
@@ -36,8 +38,9 @@ classdef Ternary_model < handle
                 T.e = params(4);
                 T.e_g0 = params(5);
                 T.u_g0 = params(6);
+                T.system_size = params(7);
                 if strcmp(ic, 'Gauss')
-                    T.x0 = params(7);
+                    T.x0 = params(8);
                 end
             end
             T.create_mesh();

@Ternary_model/create_mesh.m

@@ -3,27 +3,27 @@ function create_mesh(T)
 %   Mesh density larger for steeper gradients.
 
 % Start with left side of mesh, at r=0/x=0, with a given step size
-x = [0, 0.01];
+x = [0, 0.01/T.precision];
 
 while x<-2*T.a
     x_temp = x(end)-x(end-1);
     phi_nm1 = T.phi_tot(x(end), T.a, T.b, T.e, T.u0);
     phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0);
     if x(end) <= -T.a
-        while phi_nm1-phi_n > 0.001
+        while phi_nm1-phi_n > 0.001/T.precision
             x_temp = x_temp/2;
             phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0);
         end
     elseif x(end) > -T.a
-        while (phi_nm1-phi_n < 0.0005) && (x_temp < 0.01)
+        while (phi_nm1-phi_n < 0.0005/T.precision) && (x_temp < 0.01/T.precision)
             x_temp = 2*x_temp;
             phi_n = T.phi_tot(x(end)+x_temp, T.a, T.b, T.e, T.u0);
         end
     end
     x = [x, x(end)+x_temp];
 end
-x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 100);
-x_right = linspace(-10*T.a+0.1, T.system_size, 100);
+x_middle = linspace(-2*T.a+x(end)-x(end-1), -10*T.a, 100*T.precision);
+x_right = linspace(-10*T.a+0.1, T.system_size, 100*T.precision);
 T.x = [x, x_middle, x_right];
 end
 

@Ternary_model/plot_sim.m

@@ -2,32 +2,34 @@ function plot_sim(T, mode, nth, color)
 %PLOT_SIM Show movie of simulation
 % mode 'mov' ... show movie
 % mode 'plot' ... show plot
+si = size(T.sol);
 if nargin < 3
     nth = 1;
 end
 if strcmp(mode, 'mov')
     figure(1); hold on;
     for i = 1:nth:length(T.t)
-        cla; xlim([-T.a-1, -T.a+3]); ylim([0, 1]);
+        cla; xlim([-T.a-1, -T.a+3]); ylim([0, 2]);
         ax = gca;
         ax.FontSize = 16;
-        xlabel('position'); ylabel('volume fraction');
+        xlabel('position'); ylabel('probability [not normalized]');
         plot(T.x, T.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 2, ...
             'LineStyle', '--');
         plot(T.x, T.sol(i, :), 'LineWidth', 2);
         text(abs(-T.a+2), 0.7, num2str(T.t(i)));
         shg
         pause();
-    %     print([num2str(i),'.png'],'-dpng')
+%         print([num2str(i),'.png'],'-dpng')
     end
 elseif strcmp(mode, 'plot')
-%     figure; hold on;
-    xlim([0, 2]); ylim([0, 1]);
+%     figure;
+    hold on;
+    xlim([4, 8]); ylim([0, 2]);
     ax = gca;
     ax.FontSize = 26;
     xlabel('x [\mum]'); ylabel('volume fraction'); 
-    plot(T.x, T.sol(1:nth:150, :), 'LineWidth', 1, 'Color', color);%[0.5294, 0.8, 0.98]);
-    plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 6,...
+    plot(T.x, T.sol(1:nth:si(1), :), 'LineWidth', 1, 'Color', color);%[0.5294, 0.8, 0.98]);%
+    plot(T.x, Ternary_model.phi_tot(T.x, T.a, T.b, T.e, T.u0), 'LineWidth', 1,...
     'LineStyle', '--', 'Color', [0.97, 0.55, 0.03]);
 end
 end

@Ternary_model/solve_tern_frap.m

@@ -34,7 +34,7 @@ elseif strcmp(ic, 'Gauss')
     D_m = 1*(1-u0-e); % to make equal to ternary FRAP
     D_p = 1*(1-u0+e);
     ga = (u0-e)/(u0+e);
-    p_out = @(D_p, D_m, ga, x0, x, t) 1./(4*sqrt(D_p*pi*t))*...
+    p_out = @(D_p, D_m, ga, x0, x, t) 1./(sqrt(4*D_p*pi*t))*...
                 (exp(-(x+x0).^2./(4*D_p*t))*(ga*sqrt(D_p)-sqrt(D_m))./...
                 (ga*sqrt(D_p)+sqrt(D_m))+exp(-(x-x0).^2./(4*D_p*t)));
     if x > a

prob_laplace.m

 %% Solve Fokker Planck equation for given parameter set
-clear all; clc;
-nu = 10^-6;
-chi = 7.7/3;
-b = nu^(1/3)*sqrt(chi/(chi-2));
-e = sqrt(3/8*(chi-2));
-T = Ternary_model(0, 'Gauss', [-6, b, 0.5, e, 0.16, 0.2, 7],...
-                  linspace(0.0, 400, 5000));
-T.t = [0, 0.1, 1, 9.9];
+% Define parameters for tanh, according to Weber, Zwicker, Julicher, Lee.
+b = @(chi, nu) nu^(1/3)*sqrt(chi/(chi-2));
+e = @(chi) sqrt(3/8*(chi-2));
+%%
+T = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
+                  0.2, 300, 7], [0, 0.1, 1, 9.9], 0.5);
 T.solve_tern_frap()
 %%
+T1 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
+                   0.2, 300, 7], [0, 0.1, 1, 9.9], 1);
+T1.solve_tern_frap()
+%%
+T2 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
+                   0.2, 300, 7], [0, 0.1, 1, 9.9], 2);
+T2.solve_tern_frap()
+%%
+tic
+T3 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
+                   0.2, 300, 7], [0, 0.1, 1, 9.9], 5);
+T3.solve_tern_frap()
+toc
+%%
+T4 = Ternary_model(0, 'Gauss', [-6, b(7.7/3, 10^-6), 0.5, e(7.7/3), 0.16,...
+                   0.2, 300, 7], [0, 0.1, 1, 9.9], 10);
+T4.solve_tern_frap()
+%%
 csvwrite('x_X73_1.csv', T.x+T.a);
 csvwrite('sol_X73_1.csv', T.sol(:, :));
 %% Frank's solution to the transfer/rate problem via Laplace transform
@@ -35,7 +51,21 @@ for i = 1:length(T.t)
     shg; pause();
 end
 %% Check integral of solution, should be mass conserving and sum to 1.
-[~, ind] = min(abs(T.x-50)); % integrate up to 50, to avoid right boundary
-bins = (T.sol(:, 1:ind-1)+T.sol(:, 2:ind))/2;
-bin_size = diff(T.x(1:ind));
-sum(bins(1, :).*bin_size)
\ No newline at end of file
+[~, ind] = min(abs(T3.x-50)); % integrate up to 50, to avoid right boundary
+bins = (T3.sol(:, 1:ind-1)+T3.sol(:, 2:ind))/2;
+bin_size = diff(T3.x(1:ind));
+sum(bins(4, :).*bin_size)
+%%
+figure; hold on;
+T.plot_sim('plot', 1, 'blue');
+T1.plot_sim('plot', 1, 'magenta')
+T2.plot_sim('plot', 1, 'red');
+T3.plot_sim('plot', 1, 'green');
+T4.plot_sim('plot', 1, 'k');
+%% Check whether partitioning factor is kept throughout simulation.
+% This should work for steep boundaries.
+for i = 1:4
+pks_min = findpeaks(-T4.sol(i, :));
+pks_max = findpeaks(T4.sol(i, :));
+pks_max(1)/pks_min(1)
+end
\ No newline at end of file

ternary_frap.m

@@ -6,7 +6,6 @@ figure; % Shouldn't be symmetric, different SS for phi_u and phi_b??
 plot(min(T.u0+T.e-T.sol(:, 1:ind-1), [], 2)./max(T.u0-T.sol(:, ind+1:end), [], 2),...
     'LineWidth', 2)
 make_graph_pretty('time', 'ratio outside/inside', 'bleached component');
-%% Plot flux at boundary
 
 %% Figures
 % Time course
@@ -232,3 +231,22 @@ Ts{4}.plot_sim('plot', 10, 'yellow')
 Ts{5}.plot_sim('plot', 10, 'magenta')
 Ts{end}.plot_sim('plot', 10, 'green')
 xlim([0, 2]); ylim([0, 1]);
+%%
+nu = 10^-5;
+chi = 7/3;
+b = nu^(1/3)*sqrt(chi/(chi-2));
+e = sqrt(3/8*(chi-2));
+b = 0.05;
+e = 0.4;
+e_g0 = 0.05:0.05:0.4;
+u_g0 = 0.05:0.05:0.4;
+T = cell(length(e_g0), length(u_g0));
+tic
+parfor i = 1:length(e_g0)
+    for jj = 1:length(u_g0)
+    T(i, jj) = Ternary_model(2, 'FRAP', [-1, b, 0.5, e, e_g0(i), u_g0(jj), 40],...
+                      linspace(0.0, 100, 500));
+    T(i, jj).solve_tern_frap();
+    end
+end
+toc
\ No newline at end of file