densityVolume1D/quant_error_kd.m at main · JohnsonBiophysicsLab/densityVolume1D · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
%A0 is copies/nm.
%kf in nm/us
%kb in /us!! (not s-1!)
%Dtot in nm2/us
%sigma in nm
%time in us.
%time=logspace(-4, 7, 2000);
function[Attot, At]=quant_error_kd(A0, ka, kb, Dtot, sigma, time)


%time=logspace(-4, 7, 2000);
%solve the Smoluchoski model A+A<->A2
[Attot, res] = Rev1D_MRE_AA(A0,ka,kb,Dtot,time);

%solve the reversible A+A<->A2
kD=3*Dtot/(1/A0-sigma)
kon = (1/ka + 1/kD)^(-1)
koff=kb/ka*kon

if(kb>0)
    [At]=riccati(A0, -1, kon, koff, time);
else
    At=A0./(1+2*kon*A0*time);
end

fnum = 31;
try
    close(fnum)
    f=figure(fnum)
catch
    f=figure(fnum)
end
ax7=axes('Parent',f,'FontSize',20,'LineWidth',1,'XScale','log','YScale','linear');
hold(ax7)
plot(time/1e6, Attot(:,2)/A0,'k-','LineWidth',2)
plot(time/1e6, At/A0, 'r--','LineWidth',2)
xlabel('time (s)')
ylabel('A(t)/A(0)')
xlim([5e-6, 20])
ylim([0 1])
legend('Smoluchowski','k_{on}^{\rho}')

fnum=fnum+1
try
    close(fnum)
    f=figure(fnum)
catch
    f=figure(fnum)
end
%Aeq=Attot(end,2)
display('Aeq')
%Aeq
Aeq=Aeq_bimolecular(A0, -1, 0, kb/ka, 1/602.2)
ax7=axes('Parent',f,'FontSize',20,'LineWidth',1,'XScale','log','YScale','log');
hold(ax7)
plot(time/1e6, (Attot(:,2)-Aeq)/(A0-Aeq),'k-','LineWidth',2)
plot(time/1e6, (At-Aeq)/(A0-Aeq), 'r--','LineWidth',2)
xlabel('time (s)')
ylabel('(A(t)-Aeq)/(A(0)-Aeq)')
xlim([5e-6, 20])
ylim([0 1])
legend('Smoluchowski','k_{on}^{\rho}')
%Measure the relative error between the two curves, treating Smoluchowski
%as 'correct'.
abserr=abs(At'-Attot(:,2))/A0
relerr=abserr*A0./Attot(:,2);
[maxRelErr, relIndex]=max(relerr)
[maxAbsErr, absIndex]=max(abserr)
timeErr=time(absIndex)/1e6
timeRelErr=time(relIndex)/1e6
%define kappa

kappa = ka/kD
display('Aeq')
Aeq