fem0 matlab源码程序 - matlab其它源码

下载频道> 资源分类> matlab源码> 其它源码> fem0 matlab源码程序

标题：fem0 matlab源码程序
分享到：

所属分类: 其它源码资源类型：程序源码文件大小: 5.17 KB 上传时间: 2019-06-17 22:33:39 下载次数: 5 资源积分：1分提供者: zhangsan456 fem0

内容：
fem0 matlab源码程序

% program fem0.m

% This solves the PDE problem (D) with homogeneous Dirichlet

% conditions using piecewise linear finite elements.

clear

format long

[nodes,elements] = initial_mesh;

disp('Initial mesh is displayed in Figure Window')

visualise(nodes,elements)

disp('Type Return to Continue')

pause

nref = input('Type how many times you want to refine the initial mesh: ');

disp('Refinement Level...')

for i = 1:nref

[nodes,elements] = refine(nodes,elements);

disp(i)

end

[nnodes,m] = size(nodes);

disp('Refined mesh is displayed in Figure Window (if not too large)')

if (nnodes < 1000)

visualise(nodes,elements)

end

disp('Type Return to Continue')

pause

% First assemble the element stiffness matrices

elt_matrices = elt_stiffness(elements,nodes);

disp('Element stiffness matrices computed')

% Next assemble the global stiffness matrix

% (including all boundary nodes)

Ahat = global_stiffness(elt_matrices,elements,nodes);

disp('Global stiffness matrix computed')

N0 = find(nodes(:,3)~=1); % Finds the indices of the

% interior (non-Dirichlet) nodes

A = Ahat(N0,N0); % Select the appropriate blocks of Ahat

% to get the matrix A from lectures

% Assemble the right hand side vector (the function in the RHS of

% (D) has to be specified in source_term)

f = rhs('source_term',N0,elements,nodes);

U = A\f; % Find coefficient vector U of finite element solution

% at the degrees of freedom

% Finally for the purpose of drawing graphs of the solution

% and for computing errors, it is useful to assemble an extended

% vector Uhat which contains the elements of U at the degrees of

% freedom together with the Dirichlet boundary condition at the

% boundary nodes. The ordering of the elements of Uhat coincides

% with the ordering of the nodes.

Uhat = zeros(nnodes,1); % Define an extended vector

Uhat(N0) = U; % put in the elements of U at the degrees of freedom

% Plot the solution using the MATLAB function surf

% It is first necessary to interpolate Uhat on a uniform mesh.

% We use the MATLAB functions meshgrid and griddata to do this.

h_unif = 1/(2^(nref+1));

[XI,YI] = meshgrid([min(nodes(:,1)):h_unif:max(nodes(:,1))],...

[min(nodes(:,2)):h_unif:max(nodes(:,2))]);

ZI = griddata(nodes(:,1),nodes(:,2),Uhat,XI,YI,'linear');

hold on

surf(XI,YI,ZI);

hold off

disp('Numerical solution is displayed in Figure Window')

disp('Type Return to Continue')

pause

% Evaluate the exact solution at the nodes

x1 = nodes(:,1);

x2 = nodes(:,2);

Uex = x1.*(1-x1).*exp(x1).*x2.*(1-x2);

% Calculate error at nodes and plot it

err = abs(Uhat - Uex);

ZI = griddata(nodes(:,1),nodes(:,2),err,XI,YI,'linear');

surf(XI,YI,ZI);

disp('Error is displayed in Figure Window')

disp('Type Return to Continue')

pause

% Calculate a measure for the error (i.e. the energy norm of the difference

% between the FE solution u_h and the interpolant of the exact solution u)

anorm_err = sqrt((Uhat-Uex)'*Ahat*(Uhat-Uex))

文件列表（点击上边下载按钮，如果是垃圾文件请在下面评价差评或者投诉）:

fem0/elt_stiffness.m
fem0/fem0.m
fem0/global_stiffness.m
fem0/initial_mesh.m
fem0/mu.m
fem0/refine.m
fem0/rhs.m
fem0/source_term.m
fem0/visited.m
fem0/visualise.m

关键词： fem0 matlab源码程序

matlab源码下载排行

回到顶部

友情链接

matlab源码网站源码编程语言网络编程手机通讯编程 KUNAG 官网 KUNAG 电商更多

: 联系方式| 版权声明| 招聘信息| 广告服务| 银行汇款| 法律顾问| 兼职技术| 付款方式| 关于我们|; 网站客服程序员兼职招聘 411241577@qq.com; 沪ICP备19040327号-3; 公安备案号：沪公网安备 31011802003874号; 库纳格流体控制系统（上海）有限公司版权所有; Copyright © 1999-2014, GUSUCODE.COM, All Rights Reserved

关键词： fem0 matlab源码程序

相关推荐

matlab源码下载排行

友情链接