Mercurial > repos > public > sbplib

comparison +scheme/Schrodinger.m @ 65:33f0654a2413

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Fixed mistakes in Schrodinger scheme. BC are now working.
author Jonatan Werpers <jonatan@werpers.com>
date Wed, 18 Nov 2015 17:15:02 +0100
parents 48b6fb693025
children 446d67a49cd8
comparison
equal deleted inserted replaced
64:7067bf8adbfa 65:33f0654a2413
1 classdef Schrodinger < noname.Scheme 1 classdef Schrodinger < scheme.Scheme
2 properties 2 properties
3 m % Number of points in each direction, possibly a vector 3 m % Number of points in each direction, possibly a vector
4 h % Grid spacing 4 h % Grid spacing
5 x % Grid 5 x % Grid
6 order % Order accuracy for the approximation 6 order % Order accuracy for the approximation
8 D % non-stabalized scheme operator 8 D % non-stabalized scheme operator
9 H % Discrete norm 9 H % Discrete norm
10 M % Derivative norm 10 M % Derivative norm
11 alpha 11 alpha
12 12
13 D2
13 Hi 14 Hi
14 e_l 15 e_l
15 e_r 16 e_r
16 d1_l 17 d1_l
17 d1_r 18 d1_r
41 V_vec = V(x); 42 V_vec = V(x);
42 else 43 else
43 V_vec = x*0 + V; 44 V_vec = x*0 + V;
44 end 45 end
45 46
46 V_mat = spdiags(V,0,m,m); 47 V_mat = spdiags(V_vec,0,m,m);
47 48
48 49 obj.D = 1i * obj.D2 - 1i * V;
49 D = 1i * obj.D2 - 1i * V;
50 50
51 obj.m = m; 51 obj.m = m;
52 obj.h = h; 52 obj.h = h;
53 obj.order = order; 53 obj.order = order;
54 54
55 obj.D = alpha*obj.D2;
56 obj.x = x; 55 obj.x = x;
57 end 56 end
58 57
59 58
60 % Closure functions return the opertors applied to the own doamin to close the boundary 59 % Closure functions return the opertors applied to the own doamin to close the boundary
63 % type is a string specifying the type of boundary condition if there are several. 62 % type is a string specifying the type of boundary condition if there are several.
64 % data is a function returning the data that should be applied at the boundary. 63 % data is a function returning the data that should be applied at the boundary.
65 % neighbour_scheme is an instance of Scheme that should be interfaced to. 64 % neighbour_scheme is an instance of Scheme that should be interfaced to.
66 % neighbour_boundary is a string specifying which boundary to interface to. 65 % neighbour_boundary is a string specifying which boundary to interface to.
67 function [closure, penalty] = boundary_condition(obj,boundary,type,data) 66 function [closure, penalty] = boundary_condition(obj,boundary,type,data)
68 default_arg('type','neumann'); 67 default_arg('type','dirichlet');
69 default_arg('data',0); 68 default_arg('data',0);
70 69
71 [e,d,s] = obj.get_boundary_ops(boundary); 70 [e,d,s] = obj.get_boundary_ops(boundary);
72 71
73 switch type 72 switch type
74 % Dirichlet boundary condition 73 % Dirichlet boundary condition
75 case {'D','d','dirichlet'} 74 case {'D','d','dirichlet'}
76 tau = -1i* s * d; 75 tau = s * 1i*d;
77
78 closure = obj.Hi*tau*e'; 76 closure = obj.Hi*tau*e';
79 77
80 pp = obj.Hi*p;
81 switch class(data) 78 switch class(data)
82 case 'double' 79 case 'double'
83 penalty = pp*data; 80 penalty = -obj.Hi*tau*data;
84 case 'function_handle' 81 case 'function_handle'
85 penalty = @(t)pp*data(t); 82 penalty = @(t)-obj.Hi*tau*data(t);
86 otherwise 83 otherwise
87 error('Wierd data argument!') 84 error('Wierd data argument!')
88 end 85 end
89 86
90 % Unknown, boundary condition 87 % Unknown, boundary condition
97 % u denotes the solution in the own domain 94 % u denotes the solution in the own domain
98 % v denotes the solution in the neighbour domain 95 % v denotes the solution in the neighbour domain
99 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary); 96 [e_u,d_u,s_u] = obj.get_boundary_ops(boundary);
100 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary); 97 [e_v,d_v,s_v] = neighbour_scheme.get_boundary_ops(neighbour_boundary);
101 98
102 a = s/2 * 1i ; 99 a = s* 1/2 * 1i ;
103 b = - a'; 100 b = a';
104 101
105 tau = b*d_u; 102 tau = b*d_u;
106 sig = a*e_u; 103 sig = a*e_u;
107 104
108 closure = obj.Hi * (tau*e_u' + sig*d_u'); 105 closure = obj.Hi * (tau*e_u' + sig*d_u');