How to incorporate v and w
2 次查看(过去 30 天)
显示 更早的评论
%% u_t = D1 * u_xx, u(0,t) = u(1,t) = 0, u(x,0)= cos(pi*(x-0.5)), D1 = 0.1.
% % v_t = D2 * v_xx + a*w , v(0,t) = v(1,t) = 0, v(x,0)= sin(pi*(x-0.5)),D2 = 0.2, a = 1.5.
% % w_t = D3 * w_xx + c*u + b*v, w(0,t) = w(1,t) = 0, w(x,0)= e^(2*x),D3 = 0.3, b = 2, c = 3.
%subdivisions space /time
ht = 0.01; Tmax = 1.2; nx = 33; hx = 1/(nx-1); D1 = 0.1; x = [0:hx:1]';
%matrices
K = stiff(D1,hx,nx); M = mass(1/ht,hx,nx);A = M + K;
% disp(A(1))
%boundary conditions by deleting rows
A(nx,:) = []; A(:,nx) = []; A(1,:) = []; A(:,1) = [];
%creation
R = chol(A);
%initial condition
u = cos(pi*(x - 0.5 ));
%time step loop
k = 0;
while (k*ht < Tmax)
k = k+1;
%compute the right-hand side
b = M*u; b(nx) = []; b(1) = [];
%solve
u = zeros(nx,1); u(2:nx-1) = R\(R'\b);
figure(11)
plot(x,u(:,1),'Linewidth',1.5)
xlabel \bfx
ylabel \bfu
% % % % hold on
end
function R = stiff(nu,h,n)
[m1,m2] = size(nu);
if (length(nu) == 1)
ee = nu*ones(n-1,1)/h; e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
else
ee = 0.5*(nu(1:n-1)+nu(2:n))/h; e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
end
R = spdiags([-e1 e -e2],-1:1,n,n);
return
end
function M = mass(alpha,h,n)
[m1,m2] = size(alpha);
if(length(alpha) == 1)
ee = alpha*h*ones(n-1,1)/6;
e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
else
ee = h*(alpha(1:n-1) + alpha(2:n))/12;
e1 = [ee;0]; e2 = [0;ee]; e = e1+e2;
end
M = spdiags([e1 2*e e2],-1:1,n,n);
return
end
0 个评论
回答(0 个)
另请参阅
类别
在 Help Center 和 File Exchange 中查找有关 Workspace Variables and MAT-Files 的更多信息
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!