1

cpl11001 · Dec 14, 2017 · 5ce8415 · 5ce8415
1 parent 66a5d7d
commit 5ce8415
Showing 1 changed file with 11 additions and 7 deletions.
diff --git a/SE_diff.m b/SE_diff.m
@@ -20,10 +20,10 @@ function [pw_se,w] = SE_diff(T,P,n)
 
     % Left side of the equation
     left_side = zeros(n+1,n+1);
-    for  i = 1:n+1
-
+    for  i = 1:(n+1)^2
         left_side(i) = P * h^2*w_bar(i);
     end
+    L = sum(sum(left_side));
 
     % Compute dw_dx
     dx = zeros(n+2,n+1);
@@ -37,7 +37,7 @@ function [pw_se,w] = SE_diff(T,P,n)
         dw_dx(i,j) = (dx(i,j)+dx(i+1,j))/2;
         end 
     end 
-disp(dw_dx)
+
     % Compute dw_dy
     dy = zeros(n+1,n+2);
     for i = 2:n+2
@@ -50,10 +50,14 @@ disp(dw_dx)
         dw_dy(i,j) = (dy(i,j)+dy(i,j+1))/2;
         end 
     end 
-   disp(dw_dy)
-    pw_se = zeros(n+1,n+1);
+
+    right_side = zeros(n+1,n+1);
 
     for i = 1:(n+1)^2
-        pw_se(i) = (1000*10^3*0.3*10^-3*h^2)/(2*(1-0.31^2))*((dw_dx(i))^4/4+(dw_dy(i))^4/4+0.5*(dw_dx(i))^2+0.5*(dw_dy(i))^2);
+        right_side(i) = (0.25*(dw_dx(i))^4+0.25*(dw_dy(i))^4+0.5*(dw_dx(i)*dw_dy(i))^2);
     end 
-
+    R = sum(sum(right_side));
+    R = (1000*10^3*0.3*10^-3*h^2)/(2*(1-0.31^2))*R;
+
+    pw_se = R - L;
+end