From: Lukas Jiriste Date: Mon, 17 Jun 2024 13:10:50 +0000 (+0200) Subject: Add everything written so far X-Git-Url: https://git.ljiriste.work/?a=commitdiff_plain;h=refs%2Fheads%2Fold;p=Uchytilap.git Add everything written so far As this project was worked on before I started using versioning software, I think it best to first commit everything and then maybe fabricate an history (in alternative branch)? --- diff --git a/Kalibrace.m b/Kalibrace.m new file mode 100644 index 0000000..a5c1146 --- /dev/null +++ b/Kalibrace.m @@ -0,0 +1,189 @@ +clear all, close all +i=1; +mezvak=0.07; +cit=50; +k=50; + +Vb=36.81; +KK=616.26;%612.47 +Va=16.26; +Vv=33.87; +Vven=Va+Vb+KK; + +[fnamecelk,path]=uigetfile({'*.txt';'*.*'},'File Selector','MultiSelect','on'); +fnamecelk=cellstr(fnamecelk); + + +for i=1:length(fnamecelk) + dataV{i}=importdata([path '/' fnamecelk{i}], '\t', 2); + mn=menu(['Jaký objem je souborem ' fnamecelk{i} ' kalibrován?'],'Horní(V1)','Spodní(V2)'); + P{i,1}=dataV{i}.data(:,4-mn)'; + P{i,4}=dataV{i}.data(:,mn+1)'; +%hledání skoků + vaku=mean(P{i,1}(P{i,1}<0.01)); + P{i,4}(P{i,1}max(P{i,3})/cit)=1; + P{i,3}(P{i,3}0) + smaz=[smaz extind(l)];%zapisuje "duální" hodnoty + end + end + for l=smaz + extind(extind==l)=[]; + extind(extind>l)=extind(extind>l)-1; + P{i,1}(l)=[]; + P{i,4}(l)=[]; + P{i,2}(l)=[]; + P{i,3}(l)=[]; + end + %ještě jednou z druhé strany + smaz=[]; + for l=((1:length(extind)-1)) + if (extind(l+1)-extind(l)0) + smaz=[smaz extind(l)];%zapisuje "duální" hodnoty + end + end + for l=flip(smaz) + extind(extind==l)=[]; + extind(extind>l)=extind(extind>l)-1; + P{i,1}(l)=[]; + P{i,2}(l)=[]; + P{i,3}(l)=[]; + P{i,4}(l)=[]; + end + %a ještě jednou jen pro nárůst tlaku + smaz=[]; + minind=extind(P{i,3}(extind)==-1); + for l=flip((1:length(minind)-1)) + if minind(l+1)-minind(l)l)=extind(extind>l)-1; + P{i,1}(l)=[]; + P{i,2}(l)=[]; + P{i,3}(l)=[]; + P{i,4}(l)=[]; + end + P{i,1}(end-10:end)=[]; + P{i,2}(end-10:end)=[]; + P{i,3}(end-10:end)=[]; + P{i,4}(end-10:end)=[]; + + + figure(2) + plot(P{i,3}) + figure(1) + while extind(end)>length(P{i,3}) + extind(end)=[]; + end + if extind(1)==1 + P{i,1}(1:extind(1))=[]; + P{i,2}(1:extind(1))=[]; + P{i,3}(1:extind(1))=[]; + P{i,4}(1:extind(1))=[]; + extind=extind-extind(1); + extind(1)=[]; + end +% extind=find(abs(P{i,3})==1); + + %Ruční kontrola + kontr=1; + while kontr~=0 + plot(P{i,1},'blue') + hold on + for l=extind + if P{i,3}(l)==1 + plot([l l],[0 max(P{i,1})],'--green') + elseif P{i,3}(l)==-1 + plot([l l],[0 max(P{i,1})],'--red') + else + error(['vector extind missaligned at position ' int2str(l)]) + end + end + kontr=menu('Byla detekce bezchybná?','Ano','Ne, chci odebrat čáru','Ne, chci přebarvit čáru', 'Ne, chci přidat červenou čáru')-1; + if kontr==1 + %msgbox('Klikni na čáry označující špatně, potvrď enterem.') + [xspa,~]=ginput; + for l=xspa' + [~,ind]=min(abs(l-extind)); + extind(ind)=[]; + end + elseif kontr==2 + [xspa,~]=ginput; + for l=xspa' + [~,ind]=min(abs(l-extind)); + P{i,3}(extind(ind))=-1*P{i,3}(extind(ind)); + end + elseif kontr==3 + [xspa,~]=ginput; + for l=xspa' + + end + end + hold off + end + close all + vecdel=extind-[0 extind(1:end-1)]; + Y=mat2cell(P{i,1},1,[vecdel length(P{i,1})-extind(end)]); + + plot(P{i,1},'blue') + hold on + %třídění vektorů + extind=[0 extind length(P{i,1})]; + for j=1:length(Y) + iter=0; + X=(extind(j)+1:extind(j+1)); + Hran{i}(j,:)=[X(1) X(end)]; + suctv=[1 -1]; + while max(suctv)>2*mean(suctv) + iter=iter+1; + Koef{i}(j,:)=polyfit(X,Y{j},1); + Yodh=polyval(Koef{i}(j,:),X); + suctv=(Y{j}-Yodh).^2; + X(suctv>2*mean(suctv))=[]; + Y{j}(suctv>2*mean(suctv))=[]; + end + plot(Hran{i}(j,:),Koef{i}(j,1)*Hran{i}(j,:)+Koef{i}(j,2)) + end + hold off + Phran{i}=[]; + for j=1:length(Y)-1 + if P{i,3}(Hran{i}(j,2))==1 + p{i,1}(j)=Koef{i}(j,1)*Hran{i}(j,2)+Koef{i}(j,2); + p{i,2}(j)=Koef{i}(j+1,1)*Hran{i}(j+1,1)+Koef{i}(j+1,2); + Phran{i}=[Phran{i} [p{i,1}(j); p{i,2}(j)]]; + end + end + + Vznak(i)=mn; + V{i}=((Vv+Vven)*Phran{i}(2,:)-Vv*Phran{i}(1,:)-Vven*vaku)./(Phran{i}(1,:)-Phran{i}(2,:)); + extind(1)=[]; + extind(end)=[]; + extind(P{i,3}(extind)~=1)=[]; + P{i,5}=abs(P{i,4}-P{i,1}); + delP{i}=P{i,5}(extind); + +end +for i=1:length(V) + Vprum(i)=mean(V{i}); + printf("Výsledek z měření %s:\n\tV%i = %f mm^3\n\t", fnamecelk{i}, Vznak(i), Vprum(i)); +end +V1=mean(Vprum(Vznak==1)) +V2=mean(Vprum(Vznak==2)) +uA1=(sum((V1-[V{Vznak==1}]).^2)/(length([V{Vznak==1}])-1))^0.5 +uA2=(sum((V2-[V{Vznak==2}]).^2)/(length([V{Vznak==2}])-1))^0.5 diff --git a/allmin.m b/allmin.m new file mode 100644 index 0000000..8af01fd --- /dev/null +++ b/allmin.m @@ -0,0 +1,15 @@ +function [val, ind] = allmin(array) + dim = sum(size(array) > 1); + if dim == 0 + val = array; + ind = 1; + elseif dim == 1 + [val ,ind] = min(array); + else + ind = []; + [small_arr, single_ind] = min(array, [], dim); + [val, other_ind] = allmin(small_arr); + ind = single_ind(num2cell(other_ind){:}); + ind = [other_ind ind]; + end +end diff --git a/alternating_min.m b/alternating_min.m new file mode 100644 index 0000000..a3701db --- /dev/null +++ b/alternating_min.m @@ -0,0 +1,22 @@ +function [x, hist] = alternating_min(objfun, x_0, tol = 0.01, max_iter = 10, max_subiter = 2) + input_shape = size(x_0); + arity = prod(input_shape); + mask = zeros(input_shape); + x = x_0; + iter = 0; + rel_change = tol + 1; + while (iter < max_iter && rel_change > tol) + rel_change = 0; + iter += 1; + for i = 1:arity + mask(i) = 1; + new_fun = @(single_par) objfun(mask * single_par + ~mask .* x); + old = x(i); + x(i) = steepest_grad(new_fun, x(i), tol, max_subiter); + mask(i) = 0; + if (abs((old - x(i))/old) > rel_change) + rel_change = abs((old - x(i))/old); + end + end + end +end diff --git a/altminimize.m b/altminimize.m new file mode 100644 index 0000000..a0182ba --- /dev/null +++ b/altminimize.m @@ -0,0 +1,25 @@ +function x = altminimize(objfun, x_0, tol = 0.01, max_iter = 2) + + s = size(x_0); + if s(1) ~= 1 && s(2) ~= 1 + printf("x_0 must be given as an array of parameters!\n"); + return ; + elseif (s(2) ~= 1) + x = x_0'; + else + x = x_0; + end + + Y = zeros(size(x)); + change = (tol + 1) * ones(size(x)); + while (max(abs(change)) > tol) + for i=1:length(x) + y0 = x(i); + x(i) = 0; + Y(i) = 1; + x(i) = minimize(@(y) objfun(x + y*Y), y0, tol); + change(i) = (y0 - x(i)) / min(y0, x(i)) + Y(i) = 0; + end + end +end diff --git a/basic_model.m b/basic_model.m new file mode 100644 index 0000000..44c5895 --- /dev/null +++ b/basic_model.m @@ -0,0 +1,81 @@ +function [t, x, p] = basic_model(pars) +% Unpacking x_0 + tstep = pars(1); + xstep = pars(2); + tstart = pars(3); + tend = pars(4); + D = pars(5); + P = pars(6); + l = pars(7); + d = pars(8); + V0 = pars(9); + V1 = pars(10); + V2 = pars(11); + pmax = pars(12); + + if length(pars)<16 + all = false; + else + all = pars(16); + end + + if length(pars)<15 + T = 298; + else + T = pars(15); + end + + if length(pars)<14 + pvac = 0; + else + pvac = pars(14); + end + + if length(pars)<13 + tu = 0.5; + else + tu = pars(13); + end + +% Defining helpful constants + R = 8.314; + A = pi*d^2/4; + K0 = P*R*T*A*tstep/(V1+V0)/2/xstep; + K1 = P*R*T*A*tstep/V1/2/xstep; + K2 = P*R*T*A*tstep/V2/2/xstep; + alpha = D*tstep/xstep^2; + + m1 = tu/tstep+1; + m2 = floor((tend - tu)/tstep); + n = l/xstep+1; + m = m1+m2; + +% Initialization of resulting matrix + p = zeros(n,m); + p(1,1) = pmax; + +% Defining matrices of partial differential equations + A1 = zeros(n); + A1(2:n+1:n*n-2*n-1) = -alpha; + A1(n+2:n+1:n*n-n-1) = 1+2*alpha; + A1(2*n+2:n+1:n*n-1) = -alpha; + A1(1, 1:3) = [1+3*K0 -4*K0 K0]; + A1(end, end-2:end) = [K2 -4*K2 1+3*K2]; + + A2 = zeros(n); + A2(2:n+1:n*n-2*n-1) = -alpha; + A2(n+2:n+1:n*n-n-1) = 1+2*alpha; + A2(2*n+2:n+1:n*n-1) = -alpha; + A2(1, 1:3) = [1+3*K1 -4*K1 K1]; + A2(end, end-2:end) = [K2 -4*K2 1+3*K2]; + +% Numeric solution of Fick's law + for i = 2:m1 + p(:,i) = A1\p(:,i-1); + end + for i = m1+1:m + p(:,i) = A2\p(:,i-1); + end + t = (0:tstep:tend) + tstart; + x = 0:xstep:l; +end diff --git a/batch_sim.m b/batch_sim.m new file mode 100644 index 0000000..0a5f823 --- /dev/null +++ b/batch_sim.m @@ -0,0 +1,14 @@ +clear all +close all + +default_dir = "../Laborka/Zaloha/Matlab/"; +[fname,path] = uigetfile(default_dir, "Select 1 or more experiment files", "Multiselect", "on"); +fname = cellstr(fname); +conditions.upper_vol = input("Enter the upper volume in mm^3: ")*1e-9; % 872e-9 +conditions.lower_vol = input("Enter the lower volume in mm^3: ")*1e-9; % 600e-9 +conditions.membrane_thickness = input("Enter the membrane thickness in μm: ")*1e-6; % 350e-6 +conditions.membrane_diameter = input("Enter the effective diameter of membrane in mm: ")*1e-3; % 32e-3 +conditions.membrane_mass = input("Enter the effective mass of the membrane in g: ")*1e-3; % 0.4368e-3 +for i = 1:length(fname) + simulate_cell([path fname{i}], conditions) +end diff --git a/draw_hist.m b/draw_hist.m new file mode 100644 index 0000000..67f4475 --- /dev/null +++ b/draw_hist.m @@ -0,0 +1,16 @@ +function draw_hist(hist, dims) + if length(dims) == 1 + plot(hist(dims, :)) + return + end + points = size(hist)(2); + C = [linspace(0, 1, points); ones(1, points); ones(1, points)]'; + C = hsv2rgb(C); + if length(dims) == 2 + scatter(hist(dims(1), :), hist(dims(2), :), 10, C); + elseif length(dims) == 3 + scatter3(hist(dims(1), :), hist(dims(2), :), hist(dims(3), :), 10, C); + else + printf("The dimensionality is too high to plot!\n"); + end +end diff --git a/estimateDP.m b/estimateDP.m new file mode 100644 index 0000000..0f1362c --- /dev/null +++ b/estimateDP.m @@ -0,0 +1,27 @@ +function [D, P] = estimateDP(t, p1, p2, t_start, l, d, V0, V1, V2, T) + + R = 8.314; + A = pi * d^2 / 4; + + t = t - t_start; + p1 = p1(t>=0); + p2 = p2(t>=0); + t = t(t>=0); + + M = 20; + p1 = filtfilt(ones(M,1)/M,1,p1); + p2 = filtfilt(ones(M,1)/M,1,p2); + dp1 = p1(2:end) - p1(1:end - 1); + dp2 = p2(2:end) - p2(1:end - 1); + dt = t(2:end) - t(1:end - 1); + J1 = - V1/R/T/A * dp1/dt; + J2 = V2/R/T/A * dp2/dt; + P1 = J1 * l ./ (p1 - p2); + P2 = J2 * l ./ (p1 - p2); + P = mean([P1, P2]); + + p_eq = (p1(end) + p2(end)) / 2; + n_abs = (max(p1)*V1 + min(p2)*V2 - (V1 + V2)*p_eq)/R/T; + Sol = n_abs / (A * l) / p_eq; + D = P / Sol; +end diff --git a/expl_sol_model.m b/expl_sol_model.m new file mode 100644 index 0000000..d7c4208 --- /dev/null +++ b/expl_sol_model.m @@ -0,0 +1,86 @@ +function [t, x, p] = expl_sol_model(pars) +% Unpacking x_0 + tstep = pars(1); + xstep = pars(2); + tstart = pars(3); + tend = pars(4); + D = pars(5); + P = pars(6); + l = pars(7); + d = pars(8); + V0 = pars(9); + V1 = pars(10); + V2 = pars(11); + pmax = pars(12); + + if length(pars)<16 + all = false; + else + all = pars(16); + end + + if length(pars)<15 + T = 298; + else + T = pars(15); + end + + if length(pars)<14 + pvac = 0; + else + pvac = pars(14); + end + + if length(pars)<13 + tu = 0.5; + else + tu = pars(13); + end + +% Defining helpful constants + H = P/D; + pmax = pmax * H; # from here on p is used as a concentration + R = 8.314; + A = pi*d^2/4; + K0 = P*R*T*A*tstep/(V1+V0)/2/xstep; + K1 = P*R*T*A*tstep/V1/2/xstep; + K2 = P*R*T*A*tstep/V2/2/xstep; + alpha = D*tstep/xstep^2; + + m1 = tu/tstep+1; + m2 = floor((tend - tu)/tstep); + n = l/xstep+1; + m = m1+m2; + +% Initialization of resulting matrix + p = zeros(n,m); + p(1,1) = pmax; + +% Defining matrices of partial differential equations + A1 = zeros(n); + A1(2:n+1:n*n-2*n-1) = -alpha; + A1(n+2:n+1:n*n-n-1) = 1+2*alpha; + A1(2*n+2:n+1:n*n-1) = -alpha; + A1(1, 1:3) = [1+3*K0 -4*K0 K0]; + A1(end, end-2:end) = [K2 -4*K2 1+3*K2]; + + A2 = zeros(n); + A2(2:n+1:n*n-2*n-1) = -alpha; + A2(n+2:n+1:n*n-n-1) = 1+2*alpha; + A2(2*n+2:n+1:n*n-1) = -alpha; + A2(1, 1:3) = [1+3*K1 -4*K1 K1]; + A2(end, end-2:end) = [K2 -4*K2 1+3*K2]; + +% Numeric solution of Fick's law + for i = 2:m1 + p(:,i) = A1\p(:,i-1); + end + for i = m1+1:m + p(:,i) = A2\p(:,i-1); + end + t = (0:tstep:tend) + tstart; + x = 0:xstep:l; + +p = p/H; # This converts concentration p to equilibrium pressure p + +end diff --git a/grad.m b/grad.m new file mode 100644 index 0000000..e6a7cc2 --- /dev/null +++ b/grad.m @@ -0,0 +1,14 @@ +function g = grad(fun, args, step = 1e-6) + n = size(args); + g = zeros(n); + + argdif = zeros(n); + for i = 1:max(n) + argdif(i) = args(i) * step; + f_plus = fun(args + argdif); + f_minus = fun(args - argdif); + g(i) = (f_plus - f_minus)/(2 * argdif(i)); + argdif(i) = 0; + end +end + diff --git a/halving_int.m b/halving_int.m new file mode 100644 index 0000000..005d0b3 --- /dev/null +++ b/halving_int.m @@ -0,0 +1,17 @@ +function x = halving_int(fun, a, b, prec = 1e-10) + fa = fun(a); + fb = fun(b); + x = (a + b) / 2; + fx = fun(x); + while (abs(fx) > prec) + if (fx * fa < 0) + b = x; + fb = fx; + else + a = x; + fa = fx; + end + x = (a + b) / 2 + fx = fun(x) + end +end diff --git a/hessgrad.m b/hessgrad.m new file mode 100644 index 0000000..4f87c97 --- /dev/null +++ b/hessgrad.m @@ -0,0 +1,40 @@ +function [Hess, grad, F] = hessgrad(fun, args, step = 0.00001) + n = length(args); + grad = zeros(n, 1); + Hess = zeros(n, n); + + F = fun(args); + + argdif = zeros(n, 1); + signledif = zeros(n, 2); + for i = 1:n + argdif(i) = args(i) * step; + singledif(i, 1) = fun(args + argdif); + singledif(i, 2) = fun(args - argdif); + grad(i) = (singledif(i, 1) - singledif(i, 2))/(2 * args(i) * step); + argdif(i) = 0; + endfor + + doubledif = zeros(n, n); + for i=1:n + argdif(i) = args(i) * step; + for j=i+1:n + argdif(j) = args(j) * step; + doubledif(i, j) = fun(args + argdif); + doubledif(j, i) = fun(args - argdif); + argdif(j) = 0; + endfor + argdif(i) = 0; + endfor + + for i=1:n + for j=i:n + if j==i + Hess(i, i) = (singledif(i, 1) - 2*F + singledif(i, 2))/(args(i)*step)^2; + else + Hess(i, j) = (doubledif(i, j) + doubledif(j, i) - singledif(i, 1) - singledif(i, 2) - singledif(j, 1) - singledif(j, 2) + 2*F)/(2*args(i)*args(j)*step^2); + Hess(j, i) = Hess(i, j); + endif + endfor + endfor +end diff --git a/main.m b/main.m new file mode 100644 index 0000000..ac06a79 --- /dev/null +++ b/main.m @@ -0,0 +1,271 @@ +clear all +close all + +pkg load optim +%{ +classdef Membrane + properties + D double {mustBePositive} + P double {mustBePositive} + S double {mustBePositive} + + r double {mustBePositive} + d double {mustBePositive} + A double {mustBePositive} + l double {mustBePositive} + + material char + endproperties +endclassdef + +classdef MembraneSimulation + properties + membrane Membrane + + tstep double {mustBePositive} + xstep double {mustBePositive} + + tend double {mustBePositive} + tc double {mustBeNonnegative} = 0.5 + + pmax double {mustBePositive} + pvac double {mustBeNonnegative} = 0 + + T double {mustBePositive} = 0 + + endproperties + methods + function [t, p] = simulate(); + D = membrane.D + P = membrane.P + l = membrane.l + A = membrane.A + endfunction + endmethods +endclassdef +%} + + +[fname,path]=uigetfile; +data=importdata([path fname], '\t', 2); + +datasize=size(data.data); +t=data.data(:,1)'; +p1=data.data(:,datasize(2))'; +p2=data.data(:,datasize(2)-1)'; + +pvac = min([p1 p2]); +pmax = max(p1); + +xm=find(p1>0.5*max(p1),1,'first'); +tstart = t(xm); + +l = 350e-6; +d = 32e-3; +A = d^2/4*pi; +V0 = 616e-9; +V1 = 872e-9; +V2 = 684e-9; + + + +pars_cur = [1e-3, 3.5e-5, tstart, 100, 3e-10, 3e-13, 3.5e-4, 32e-3, V0, V1, V2, pmax, 0.3]; + +function [tmod, p] = cellmod(pars) + % cellmod(tstep, xstep, tstart, tend, D, P, l, d, V0, V1, V2, pmax, tc=0.5, pvac=0, T = 298, all) + + tstep = pars(1); + xstep = pars(2); + tstart = pars(3); + tend = pars(4); + D = pars(5); + P = pars(6); + l = pars(7); + d = pars(8); + V0 = pars(9); + V1 = pars(10); + V2 = pars(11); + pmax = pars(12); + + + if length(pars)<16 + all = false; + else + all = pars(16); + end + + if length(pars)<15 + T = 298; + else + T = pars(15); + end + + if length(pars)<14 + pvac = 0; + else + pvac = pars(14); + end + + if length(pars)<13 + tc = 0.5; + else + tc = pars(13); + end + + A = d^2/4*pi; + + R = 8.314; + + m = floor(l/xstep)+1; + xstep = l/m; + + n1 = floor(tc/tstep)+1; + tstep1 = tc/n1; + + if tend>tc + n2 = floor( (tend-tc)/tstep); + tstep2 = (tend-tc)/n2; + p = zeros(n1+n2-1, m); + else + tstep2 = tstep; + p = zeros(n1, m); + end + + % criterium for numerical stability + crit = D*max([tstep1, tstep2])/xstep^2; + while crit>1 % Meaning error grows + if tstep1 > tstep2 + tstep1 = tstep1/2; + n1 = n1*2; + else + tstep2 = tstep2/2; + if tend > tc + n2 = n2*2; + end + endif + crit = D*max([tstep1, tstep2])/xstep^2; + end + + tstep1; + tstep2; + + tmod(1) = 0; + + p(1, :) = pvac; + p(1,1) = pmax; + + for j=1:n1 + p(j+1, 2:end-1) = p(j, 2:end-1) + ... + D*tstep1/xstep^2*( p(j, 1:end-2) - 2*p(j, 2:end-1) + p(j, 3:end)); + p(j+1, 1) = p(j, 1) + tstep1/2/xstep*P*R*T*A/(V0+V1)*... + ( -3*p(j, 1) + 4*p(j, 2) - p(j, 3) ); + p(j+1, m) = p(j, m) - tstep1/2/xstep*P*R*T*A/V2*... + ( p(j, m-2) - 4*p(j, m-1) + 3*p(j, m) ); + + tmod(j+1) = tmod(j) + tstep1; + end + if tend > tc + for j=n1+1:n1+n2-1 + p(j+1, 2:end-1) = p(j, 2:end-1) + ... + D*tstep1/xstep^2*( p(j, 1:end-2) - 2*p(j, 2:end-1) + p(j, 3:end)); + p(j+1, 1) = p(j, 1) + tstep2/2/xstep*P*R*T*A/V1*... + ( -3*p(j, 1) + 4*p(j, 2) - p(j, 3) ); + p(j+1, m) = p(j, m) - tstep2/2/xstep*P*R*T*A/V2*... + ( p(j, m-2) - 4*p(j, m-1) + 3*p(j, m) ); + + tmod(j+1) = tmod(j) + tstep2; + end + end + + tmod = [0 tstart tmod+tstart]; + p = [pvac*ones(2,size(p,2)); p]; + + if all<0.5 + p = [p(:,1) p(:,end)]; %Erasing all but the boundary (cell) pressures + endif +end + +function Obj = objfun(model, pars, exp) + xexp = exp(:,1); + yexp = exp(:,2:end); + + [xmod, ymod] = model(pars); + + ymod_interpol = interp1(xmod, ymod, xexp); + + valid = ~isnan(ymod_interpol) & (yexp>0); + + Obj = 1/sum(sum(valid)) * sum(sum( (yexp(valid) - ymod_interpol(valid)).^2 ./ yexp(valid) )); +end + +function grad = grad(field, pars, ind, step = 0.001) + grad = zeros(size(pars)); + temp_pars = pars; + for i = ind + j = j+1; + + temp_pars_pos = pars; + temp_pars_pos(i) = temp_pars_pos(i)*(1+step); + + val_pos = field(temp_pars_pos); + + temp_pars_neg = pars; + temp_pars_neg(i) = temp_pars_neg(i)*(1-step); + + val_neg = field(temp_pars_neg); + + grad(i) = (val_pos - val_neg)/(temp_pars_pos(i) - temp_pars_neg(i)); + endfor +end + +tolobj = 1e-4; + +curobjfun = @(pars) objfun(@cellmod, pars, [t; p1; p2]'); + + +pars_prev = pars_cur; +pars_prev(5:6) = pars_prev(5:6)*1.01; +gradF_prev = grad(curobjfun, pars_prev, [5 6]); +F_prev = curobjfun(pars_prev); + +krit = tolobj+1; + +optstep = 1; + +while F_prev>tolobj + + + gradF_cur = grad(curobjfun, pars_cur, [5 6]) + + pars_cur = pars_cur - optstep * numhessian(curobjfun, pars_cur, [5 6]) \ gradF_cur + F_cur = curobjfun(pars_cur) +%{ + while F_cur > F_prev + optstep = optstep/2; + pars_cur = pars_prev - optstep * gradF; + F_cur = curobjfun(pars_cur); + endwhile + + krit = abs(F_prev - F_cur) + F_cur + optstep + pars_cur + gradF_cur +%} + F_prev = F_cur; + gradF_prev = gradF_cur; + pars_prev = pars_cur; + +end + +%{ +p1mod = p(:,1); +p2mod = p(:,end); + +t = (0:0.01:tend); +x = (0: 10e-6 : 350e-6); + +[X,T] = meshgrid(x,t); + +surf(X,T,p) +%} diff --git a/membrane.m b/membrane.m new file mode 100644 index 0000000..7018a84 --- /dev/null +++ b/membrane.m @@ -0,0 +1,5 @@ +classdef membrane + properties + name + radius + diff --git a/objfun.m b/objfun.m new file mode 100644 index 0000000..12a83d1 --- /dev/null +++ b/objfun.m @@ -0,0 +1,20 @@ +function [E, t_calc, x, p] = objfun(model, pars, experimental) +% experimental = [t_exp, p1_exp, p2_exp] + t_exp = experimental(1, :); + p1_exp = experimental(2, :); + p2_exp = experimental(3, :); + + [t_calc, x, p] = model(pars); + + first_ind = find(t_exp >= min(t_calc), 1, 'first'); + last_ind = find (t_exp <= max(t_calc), 1, 'last'); + + t_exp = t_exp(first_ind:last_ind); + p1_exp = p1_exp(first_ind:last_ind); + p2_exp = p2_exp(first_ind:last_ind); + + p1_calc = interp1(t_calc, p(1, :), t_exp); + p2_calc = interp1(t_calc, p(end, :), t_exp); + + E = (sumsq((p1_calc - p1_exp)) + sumsq((p2_calc - p2_exp)))/(length(p1_calc) + length(p2_calc)); +end diff --git a/octave-workspace b/octave-workspace new file mode 100644 index 0000000..84c5b3b Binary files /dev/null and b/octave-workspace differ diff --git a/simulate_cell.m b/simulate_cell.m new file mode 100644 index 0000000..ab3993e --- /dev/null +++ b/simulate_cell.m @@ -0,0 +1,149 @@ +function simulate_cell(full_path, conditions) + %{ + everything is in SI base units + conditions + mandatory: + upper_vol % V1 + lower_vol % V2 + membrane_mass % m + membrane_thickness % l + membrane_diameter % diameter of the part of membrane + % that is exposed inside the aparatus + voluntary: + temperature % temp + outer_vol % volume from which the gas expands inside + %} + + %{ + figure + surf(t, x, p); + shading flat + %} + + pkg load signal + + data=importdata(full_path, '\t', 2); + datasize=size(data.data); + t_exp=data.data(:,1)'; + p1=data.data(:,datasize(2))'*1e5; + p2=data.data(:,datasize(2)-1)'*1e5; + + a=0; + b=0; + + pvac = min([p1 p2]); + pmax = max(p1); + + xm=find(p1>0.5*max(p1),1,'first'); + t_start = t_exp(xm); + + R = 8.314; + % T = 295; + if (isfield(conditions, "temperature")) + T = conditions.temperature; + else + T = 295; + endif + if (isfield(conditions, "outer_vol")) + V0 = conditions.outer_vol; + else + V0 = inf; + endif + % V0 = 616e-6; + % V1 = 872e-9; + % V2 = 600e-9; + % l = 350e-6; + % d = 32e-3; + % m = 0.4368e-3; + V1 = conditions.upper_vol; + V2 = conditions.lower_vol; + l = conditions.membrane_thickness; + d = conditions.membrane_diameter; + m = conditions.membrane_mass; + tu = 0.5; + + A = pi * d^2 / 4; + + tc = 100; + + %[D, P] = estimateDP(t_exp, p1, p2, t_start, l, d, V0, V1, V2, T); + D = 1.97093e-9; + P = 1.04363e-12; + H = P/D; + + %printf("First estimate of membrane characteristics:\tD = %g\tP = %g\n\n", D, P); + + tstep = tu/100; + xstep = l/50; + + curobjfun = @(pars) objfun(@basic_model, [tstep, xstep, t_start, 300, pars(1), pars(2), l, d, V0, V1, V2, pmax], [t_exp; p1; p2]); + + x = [D, P]; + %x = steepest_grad(curobjfun, x, 0.01, 20); + x = alternating_min(curobjfun, x, 0.01, 10); + D = x(1); + P = x(2); + printf("Final estimate of membrane characteristics:\tD = %g\tP = %g\n\n", D, P); + + [E, t, l, p] = curobjfun(x); + + printf("The error of the final estimate is %g Pa^2.\n\n", E); + + %mkdir(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes]) + + figure() + plot(t, p(1,:), 'r', t, p(end,:), 'r', t_exp, p1, 'b', t_exp, p2, 'b'); + legend('p\_calc', '', 'p\_exp', '') + ylabel('p', 'Rotation', 0) + xlabel('t') + + %print(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes '/tlak.png'],'-dpng','-r360')%,'-painters') + + J1_mem = - P * (-3*p(1,:) + 4*p(2,:) - p(3,:)) / (2*xstep); + J2_mem = - P * (3*p(end,:) - 4*p(end-1,:) + p(end-2,:)) / (2*xstep); + figure() + plot(t, J1_mem, t, J2_mem) + hold on + + before_closed = (t - t_start) < tu; + V_up = (V0 + V1) * before_closed + V1 * (~before_closed); + + J1_der = -V_up(2:end-1).*(1/A/R/T * (p(1, 3:end) - p(1, 1:end-2)) / (2 * tstep)); + J1_der(tu/tstep) = (J1_der(tu/tstep - 1) + J1_der(tu/tstep + 1)) / 2; + J1_der = [ (-V_up(1) * 1/A/R/T * (-3*p(1, 1) + 4*p(1, 2) - p(1, 3)) / (2 * tstep)) J1_der]; + J1_der = [ J1_der (-V_up(end) * 1/A/R/T * (3*p(1, end) - 4*p(1, end - 1) + p(1, end - 2)) / (2 * tstep))]; + J2_der = V2/A/R/T * (p(end, 3:end) - p(end, 1:end-2)) / (2 * tstep); + J2_der = [ (V2 * 1/A/R/T * (-3*p(end, 1) + 4*p(end, 2) - p(end, 3)) / (2 * tstep)) J2_der]; + J2_der = [ J2_der (V2 * 1/A/R/T * (3*p(end, end) - 4*p(end, end - 1) + p(end, end - 2)) / (2 * tstep))]; + plot(t, J1_der, t, J2_der) + legend('J1_x', 'J2_x', 'J1_t', 'J2_t') + ylabel('J', 'Rotation', 0) + xlabel('t') + hold off + + %print(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes '/tok.png'],'-dpng','-r360')%,'-painters') + + n = zeros(size(J1_mem)); + n(1) = tstep * (J1_mem(1) - J2_mem(1)); + for i=2:length(n) + n(i) = n(i-1) + tstep * (J1_mem(i) - J2_mem(i)); + end + n = n * A / m; + n_der = zeros(size(J1_der)); + n_der(1) = tstep * (J1_der(1) - J2_der(1)); + for i=2:length(n_der) + n_der(i) = n_der(i-1) + tstep * (J1_der(i) - J2_der(i)); + end + n_der = n_der * A / m; + n_bil = 1 / R / T * (max(p1) * V1 - p1 * V1 - p2 * V2) / m; + figure() + plot(t, n, t, n_der, t_exp, n_bil) + legend('n_x', 'n_t', 'n_{bil}') + ylabel('n/m', 'Rotation', 0) + xlabel('t') + + %print(['Obrázky/PNG/NumFick/' Membrana '_' plyn '_p0=' tlak '_' dnes '/sorpce.png'],'-dpng','-r360')%,'-painters') + + p_end = (p(1,end)*V1 + p(end,end)*V2)/(V1 + V2)/100000; + printf("Solubility of the gas in mol/kg/bar by method is:\n\tSpatial\t\tTemporal\tBilance\n\t%d\t%d\t%d\n", n(end)/p_end, n_der(end)/p_end, n_bil(end)/p_end) +end diff --git a/steepest_grad.m b/steepest_grad.m new file mode 100644 index 0000000..954e78d --- /dev/null +++ b/steepest_grad.m @@ -0,0 +1,56 @@ +function [x, hist] = steepest_grad(objfun, x_0, tol = 0.01, max_iter = 5) + s = size(x_0); + if s(1) ~= 1 && s(2) ~= 1 + printf("x_0 must be given as an array of parameters!\n"); + return ; + elseif (s(2) ~= 1) + x = x_0'; + else + x = x_0; + end + + F = objfun(x); + hist = [x; F]; + reldel = tol + 1; + iter = 0; + g = - grad(objfun, x); + alpha = 0.9 * min(abs(x./g)); + + while (abs(reldel) > tol && iter < max_iter) + iter = iter + 1; + %printf("Iteration number:\t%i\n", iter) + % Armijo-Goldstein condition + tau = 0.5; + c = 0.1; + alpha = 0.9 * min(abs(x./g)); + m = g*g'; + t = - c * m; + F_new = objfun(x + alpha * g); + while (F - F_new < alpha * t) + %printf("Decreasing alpha.\n") + alpha = alpha * tau; + if (min(abs((alpha * g) ./ x)) < sqrt(eps)) + printf("Alpha too small, reached a very sharp \"valey\".\n") + return ; + end + F_new = objfun(x + alpha * g); + end + reldel = 1 - F_new / F; + F = F_new; + hist = [hist [(x + alpha * g); F_new]]; + if (any(isnan(x + alpha * g))) + printf("NaN encountered, ending minimization!\n"); + break + end + x = x + alpha * g; + %{ + printf("\nImprovement:\t\t%f %%\n", reldel * 100); + printf("Current parameters:"); + for i = 1:length(x)-1 + printf("\t%g, ", x(i)); + end + printf("\t%g\n\n", x(end)); + %} + g = - grad(objfun, x); + end +end