From 42bf60f2cd0858a2dfbb24d053de633424e69459 Mon Sep 17 00:00:00 2001
From: Lars Hubatsch <hubatsch@pks.mpg.de>
Date: Tue, 19 Nov 2019 14:36:06 +0100
Subject: [PATCH] Analytical derivates work, parameter regime doesn't look
completely off.
---
ternary_frap.m | 46 +++++++++++++++++++++++++++++++---------------
1 file changed, 31 insertions(+), 15 deletions(-)
diff --git a/ternary_frap.m b/ternary_frap.m
index 3e7b8ac..70df44d 100644
--- a/ternary_frap.m
+++ b/ternary_frap.m
@@ -4,8 +4,8 @@
% conditions are introduced. Integration of model via pdepe.
%% Solve pde
-x = linspace(0, 40, 500);
-t = linspace(0, 5000, 50);
+x = linspace(0,500, 6000);
+t = linspace(0, 500, 50);
% opt = odeset('RelTol',1e-14, 'AbsTol', 1e-16,'MaxStep',1e-2);
% sol = pdepe(0,@flory_pde, @flory_ic, @flory_bc,x,t, opt);
sol = pdepe(0,@flory_hugg_pde, @flory_ic, @flory_bc,x,t);
@@ -13,28 +13,36 @@ sol = pdepe(0,@flory_hugg_pde, @flory_ic, @flory_bc,x,t);
%% Plotting
figure(1); hold on;
for i = 1:50
- cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.01);
+ cla; plot(x, pt(x)); plot(x, sol(i, :)); pause(0.04);
end
%% Plot and check derivatives of pt
figure; hold on;
+x = linspace(140, 170, 500);
plot(x, pt(x)); plot(x, gra_pt(x, 0.001)); plot(x, lap_pt(x, 0.001));
plot(x, gralap_pt(x, 0.001)); plot(x, laplap_pt(x, 0.001));
-
+plot(x, gra_a(x, -155, 0.5));
+plot(x, lap_a(x, -155, 0.5));
+% pt(0)
+% gra_pt(0, 0.00001)
+% lap_pt(0, 0.00001)
+% lap_a(0, -155, 0.5)
+% gralap_pt(0, 0.00001)
+% laplap_pt(0, 0.00001)
%% Function definitions for pde solver
function [c, f ,s] = flory_hugg_pde(x, t, u, dudx)
% Solve with full ternary model.
-X = 2.2;
-K = 1;
+X = 3.2;
+K = 1.5;
c = 1/(pt(x)+1);
f = dudx;
-s = u*(lap_pt(x, 0.001)-2*X*(lap_pt(x, 0.001)-...
- gra_pt(x, 0.001)^2-pt(x)*lap_pt(x, 0.001))+...
- K*(-laplap_pt(x, 0.001)+...
- gra_pt(x, 0.001)*gralap_pt(x, 0.001)+...
- pt(x)*laplap_pt(x, 0.001)))/(pt(x)+1)+...
- dudx*(-2*X*(gra_pt(x, 0.001)-pt(x)*gra_pt(x, 0.001))-...
- K*gralap_pt(x, 0.001)+K*pt(x)*gralap_pt(x, 0.001))/...
+s = u*(lap_pt(x, 0.0001)-2*X*(lap_pt(x, 0.0001)-...
+ gra_pt(x, 0.0001)^2-pt(x)*lap_pt(x, 0.0001))+...
+ K*(-laplap_pt(x, 0.0001)+...
+ gra_pt(x, 0.0001)*gralap_pt(x, 0.0001)+...
+ pt(x)*laplap_pt(x, 0.0001)))/(pt(x)+1)+...
+ dudx*(-2*X*(gra_pt(x, 0.0001)-pt(x)*gra_pt(x, 0.0001))-...
+ K*gralap_pt(x, 0.0001)+K*pt(x)*gralap_pt(x, 0.0001))/...
(pt(x)+1);
end
@@ -42,7 +50,7 @@ function u0 = flory_ic(x)
if x<5
u0 = 0.0;
else
- u0 = 0.1;
+ u0 = 0.5;
end
% u0 = pt(x)*0.5;
end
@@ -64,7 +72,7 @@ s = -u*lap_pt(x)/pt(x);
end
function p = pt(x)
- p = (tanh(-(x-5)*2)+1)/2+0.1;
+ p = (tanh(-(x-155)*2)+1)/2+0.3;
end
function gpt = gra_pt(x, delta)
@@ -84,4 +92,12 @@ end
function llpt = laplap_pt(x, delta)
llpt = (gralap_pt(x+delta, delta)-...
gralap_pt(x-delta, delta))/(2*delta);
+end
+
+function grad = gra_a(x, a, b)
+ grad = -(1-tanh(-(x+a)/b).^2)*1/b*0.5;
+end
+
+function lap = lap_a(x, a, b)
+ lap = -2*tanh(-(x+a)/b).*(1-tanh(-(x+a)/b).^2)*1/b^2*0.5;
end
\ No newline at end of file
--
GitLab