matlab有限元分析实例
MATLAB有限元分析与应用精选全文完整版
%SpringElementForces This function returns the element nodal force
%
vector given the element stiffness matrix k
%
and the element nodal displacement vector u.
2019/11/28
§2-1 弹簧元
u1=U(1:2); f1=SpringElementForces(k1,u1);
f1 = -15.0000 15.0000
u2=U(2:3); f2=SpringElementForces(k2,u2);
f2 = -15.0000 15.0000
12
§3-1 弹簧元
%
modulus of elasticity E, cross-sectional
%
area A, and length L. The size of the
%
element stiffness matrix is 2 x 2.
y = [E*A/L -E*A/L ; -E*A/L E*A/L];
2019/11/28
3.1 单元刚度矩阵的形成
function y = SpringElementStiffness(k)
%SpringElementStiffness This function returns the element stiffness %matrix for a spring with stiffness k. %The size of the element stiffness matrix is 2 x 2.
matlab有限元三角形单元编程
matlab有限元三角形单元编程
在MATLAB中进行有限元分析，可以使用其自带的有限元分析工具箱（FEATool）进行编程。
以下是一个简单的例子，演示如何使用三角形单元进行有限元分析：
1. 打开MATLAB，进入FEATool环境。
2. 创建新的有限元模型。
选择“Model”菜单下的“New Model”选项，进入“Model Builder”界面。
3. 在“Model Builder”界面中，选择“2D Triangle”单元类型，并在绘图区域中绘制出三角形网格。
4. 在“Model Builder”界面中，设置材料属性、边界条件和载荷等参数。
5. 运行有限元分析。
选择“Model”菜单下的“Solve”选项，进行有限元求解。
6. 查看结果。
选择“Model”菜单下的“Results”选项，可以查看位移、应力、应变等结果。
以上是一个简单的例子，演示了如何使用三角形单元进行有限元分析。
在实际应用中，还需要根据具体问题进行详细的建模和计算。
Matlab有限元分析操作基础
Matlab有限元分析操作基础Matlab 有限元分析20140226为了用Matlab 进行有限元分析，首先要学会Matlab 基本操作，还要学会使用Matlab 进行有限元分析的基本操作。
1. 复习：上节课分析了弹簧系统x 推导了系统刚度矩阵1122121200k k k k k k k k ----+??2. Matlab有限元分析的基本操作（1）单元划分（选择何种单元，分成多少个单元，标号）（2）构造单元刚度矩阵（列出…）（3）组装系统刚度矩阵（集成整体刚度矩阵）（4）引入边界条件（消除冗余方程）（5）解方程（6）后处理（扩展计算）3. Matlab有限元分析实战【实例1】分析：步骤一：单元划分步骤二：构造单元刚度矩阵>>k1=SpringElementStiffness(100) >>…?步骤三：构造系统刚度矩阵a) 分析SpringAssemble库函数function y = SpringAssemble(K,k,i,j)% This function assembles the element stiffness% matrix k of the spring with nodes i and j into the % global stiffness matrix K.% function returns the global stiffness matrix K% after the element stiffness matrix k is assembled. K(i,i) = K(i,i) + k(1,1);K(i,j) = K(i,j) + k(1,2);K(j,i) = K(j,i) + k(2,1);K(j,j) = K(j,j) + k(2,2);y = K;b) K是多大矩阵？今天的系统刚度矩阵是什么？因为11221212k kk kk k k k----+所以10001000200200 100200300----c) K=SpringAssemble(K,k1,1,2) function y = SpringAssemble(K,k,i,j) K(i,i) = K(i,i) + k(1,1);K(i,j) = K(i,j) + k(1,2);K(j,i) = K(j,i) + k(2,1);K(j,j) = K(j,j) + k(2,2); 1100100100100k -??=??-??10010001001000000K -??=-??K=SpringAssemble(K,k2,2,3) 1200200200200k -??10010001003002000200200K -??=---??1001000100010010030020002002000200200100200300----≠----步骤四：引入边界条件，消除冗余方程>>k=K(2:3,2:3)%构造不含冗余的方程>>f=[0;15]%构造外力列阵步骤五：解方程引例：已知1212u 31u u u +=??-=?，求 12u u 和解：类似求解KU=F ，输入下列Matlab 命令：>> K=[1 1;1,-1]>> F=[3;1]>> U=inv(K)*F>> U=K \F（继续弹簧系统求解）>>u=k \f %使用高斯消去法求解>>U=[0 ; u]%构造原方程组>>F=K*U %求出所有外力，含多余计算步骤六：后处理、扩展计算>>u1=[0;U(2)]%构造单元位移>>f1=SpringElementForces(k1,u1)%求单元1内力>>u2=[U(2) ; U(3)]%构造单元2位移>>f2=SpringElementForces(k2,u2)%求单元2内力4. 总结clck1=SpringElementStiffness(100)%创建单元刚度矩阵 1 k2=SpringElementStiffness(200)%创建单元刚度矩阵 2 K=zeros(3,3)%创建空白整体刚度矩阵K=SpringAssemble(K,k1,1,2)%按节点装入单元矩阵 1 K=SpringAssemble(K,k2,2,3)%按节点装入单元矩阵2 k=K(2:3,2:3)%构造不含冗余的方程f=[0;15]%构造外力列阵u=k\f%使用高斯消去法求解U=[0 ; u]%构造系统节点位移列阵F=K*U%求出所有外力，含多余计算u1=[0;U(2)]%构造单元位移f1=SpringElementForces(k1,u1)%求单元1内力u2=[U(2) ; U(3)]%构造单元2位移f2=SpringElementForces(k2,u2)%求单元2内力5. 练习1 Danyi 132 dan 34 3dan 35 4dan 35 dan5 54 dan6 42。
悬臂梁MATLAB有限元算例注释
用有限元法对悬臂梁分析的算例算例：如下图所示的悬臂梁，受均布载荷q ＝1N ／mm 2作用。
E ＝2．1×105N ／mm 2, μ＝0.3厚度h ＝10mm 。
现用有限元法分析其位移及应力。
梁可视为平面应力状态，先按图示尺寸划分为均匀的三角形网格，共有8×10＝80个单元，5×ll ＝55个节点，坐标轴以及单元与节点的编号如图。
将均布载荷分配到各相应节点上，把有约束的节点5l 、52、53、54、55视作固定铰链，建立如图所示的离散化计算模型。
程序计算框图：（续左）（接右）开 始 输入材料参数 计算具有代表性的单元刚阵 K<=0 将各单元刚阵按整体编号集成到整体刚阵 处理根部约束，修改【K 】【Q 】 求解[K][δ]＝[Q] 整理[δ] 并画图计算单元应力，并输出结束程序中的函数功能介绍及源代码1. LinearTriangleElementStiffness(E,NU,t,xi,yi,xj,yj,xm,ym)――该函数用于计算平面应力情况下弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)、第二个节点坐标为(xj,yj)、第三个节点坐标为(xm，ym)时的线性三角形元的单元刚度矩阵.该函数返回6×6的单位刚度矩阵k.2. LinearTriangleAssemble(K,k,i,j,m)――该函数将连接节点i,j,m的线性三角形元的单元刚度矩阵k集成到整体刚度矩阵K。
每集成一个单元，该函数都将返回2N×2N的整体刚度矩阵K.3. LinearTriangleElementStresses(E,NU,t,xi,yi,xj,yj,xm,ym,u)-- 该函数计算在平面应力情况下弹性模量为E、泊松比为NU、厚度为t、第一个节点坐标为(xi,yi)第二个节点坐标为(xj,yj)、第三个节点坐标为(xm，ym)以及单元位移矢量为u时的单元应力。
运用MATLAB对桁架单元进行有限元分析
np=3; ne=2; nload=1; nb=4; nu=0;% np为节点数，ne为单元数，nload为外力数，nb为约束数，nu 为泊松比np2=np*2; np3=np2-nb;% np2为不受约束时自由度是节点的两倍 ，np3为不受约束的节点自由度个数K=sym(zeros(np2,np2));% 定义受整体刚度空数组F=sym(zeros(np2,1));% 定义节点外力空数组KK=sym(zeros(np3,np3));% 预置自由度刚度空数组FF=sym(zeros(np3,1));% 预置自由度外力空数组syms h A P E L% 定义未知正常数为变量xy=[0,h;sqrt(3)*h,0;0,0];% 节点横纵坐标数组AA=[A;A];% 单元杆件面积load=[2,2,-P];% 荷载数组bound=[1,1,0;1,2,0;3,1,0;3,2,0];% 边界条件数组IJ=[1,2;3,2];% 各杆单元节点编号DW=zeros(1,4);for ie =1:neip=IJ(ie,1);jp=IJ(ie,2);DW(1)=ip*2-1; DW(2)=ip*2; DW(3)=jp*2-1; DW(4)=jp*2; % 给单元节点横纵方向编号x1=xy(ip,1); x2=xy(jp,1); y1=xy(ip,2); y2=xy(jp,2); % 杆单元节点x，y坐标L=sqrt((x2-x1)^2+(y2-y1)^2);% 杆单元长度c=(x2-x1)/L; s=(y2-y1)/L;% c为余弦， s为正弦T=[c,s,0,0;-s,c,0,0;0,0,c,s;0,0,-s,c];% 坐标转换矩阵（将局部坐标转换为整体坐标）A1=AA(ie);ke=[E*A1/L,0,-E*A1/L,0;0,0,0,0;-E*A1/L,0,E*A1/L,0;0,0,0,0];k=T.'*ke*T;% 将局部坐标下单元刚度转换为整体坐标下单元刚度 (转置后面加.可以去掉结果中的虚数)for i =1:4i1=DW(i);for j =1:4j1=DW(j);K(i1,j1)=K(i1,j1)+k(i,j);% 集成整体刚度矩阵endendendfor i =1:nloadi1=(load(i,1)-1)*2+load(i,2);F(i1,1)=F(i1,1)+load(i,3);% 将荷载按节点方向代码累加，计入外力荷载endFuu=sym(zeros(np2,1));% 节点位移NR=zeros(np2,1);for i =1:nbi1=(bound(i,1)-1)*2+bound(i,2);NR(i1)=-i1;% 将有约束的节点的横纵节点编号挑出来变为负数uu(i1)=bound(i,3);% 有约束的四个方向位移为0endj=0;for i =1:np2i1=NR(i);if i1==0j=j+1;NR(i)=j;endend% 此循环目的是为了将没有约束的节点以及方向挑出来for i =1:np2i1=NR(i);if i1>0FF(i1)=F(i);end% 将没有约束的节点方向的外力储存下来for j =1:np2j1=NR(j);if (i1>0 & j1>0)KK(i1,j1)=K(i,j);end% 将没有约束的节点的刚度储存下来endendKKFFfor i =1:np2i1=NR(i);for j =1:np2j1=NR(j);if (i1>0 & j1<0)jj1=abs(j1);FF(i1)=FF(i1)-K(i,jj1)*uu(jj1);end% 由于有约束方向的位移为0，所以无位移的地方FF[i1]不变endendu=sym(zeros(np2,1));u=KK\FF% 求解线性方程for i =1:np2i1=NR(i);if i1>0uu(i)=u(i1);end% 将求解的无约束方向的位移放进总的位移数组相应的位置中去enduudisp('位移')for ip =1:npstr=[ip;uu(ip*2-1);uu(ip*2)];disp(str)% 输出节点位移endfor ie =1:neip=IJ(ie,1);jp=IJ(ie,2);DW(1)=ip*2-1; DW(2)=ip*2; DW(3)=jp*2-1; DW(4)=jp*2; x1=xy(ip,1); x2=xy(jp,1); y1=xy(ip,2); y2=xy(jp,2);L=sqrt((x2-x1)^2+(y2-y1)^2);c=(x2-x1)/L; s=(y2-y1)/L;T=[c,s,0,0;-s,c,0,0;0,0,c,s;0,0,-s,c];A1=AA(ie);ke=[E*A1/L,0,-E*A1/L,0;0,0,0,0;-E*A1/L,0,E*A1/L,0; 0,0,0,0];d=sym(zeros(4,1));for i=1:4d(i)=uu(i);endde=T*d;% 将整体坐标位移转换为局部坐标位移Fe=ke*de;% 杆单元力for j=1:4N=Fe((ie-1)*2+1);% 杆单元轴力sigma=N/A1;enddisp('单元，轴力，应力')str=[ie;N;sigma];disp(str)endF3=K*uu;R=F3-F;disp('支反力');for i =1:nbi1=(bound(i,1)-1)*2+bound(i,2);str=[bound(i,1);bound(i,3);simplify(R(i1))]; disp(str)end。
matlab编程平面梁问题有限元分析
前处理程序clear;clcXY = [1, 0, 02, 2, 03, 2, -2]; %XY为N行3列，N为节点总数ELB = [1, 1, 2, 12, 2, 3, 1]; %ELB为MB行4列，MB为单元总数b = 4.5e-2;h = 2e-1; %截面宽b和高hA1 = b*h;IZ1 = 3e-5; %惯性矩EAIZ = [200e9, A1, IZ1]; %弹性模量与截面积和惯性矩ELPQ = [1, 1, -2e4, 22, 1, -1e4, 2]; %ELPQ第一列为受载荷单元号，第二列为长度C（见下表），三列为载荷大小，四列为载荷类型SU = [1,02, 03,0]; %SU第一列为已知节点位移对应的方程号，第2列为节点位移数值，约束即等于0子程序1function [PO, ii, jj] = dxjd(ELB, XY, ELPQ1)PO = zeros(6,1);k = ELPQ1(1);jj = ELB(k, 3);dxy = [XY(ii, 2), XY(ii, 3)XY(jj, 2), XY(jj, 3)];DY = dxy(2,2) - dxy(1,2);DX = dxy(2,1) - dxy(1,1);L = sqrt(DX^2 + DY^2);C = ELPQ1(2);Q = ELPQ1(3);type = ELPQ1(4);switch typecase 1PO(2) = PO(2) + 0.5*Q*C*(2-2*C^2/L^2 + C^3/L^3);PO(3) = PO(3) + Q*C^2*(6-8*C/L+3*C^2/L^2)/12;PO(5) = PO(5)+Q*C-PO(2);PO(6) = PO(6) - Q*C^3*(4-3*C/L)/12/L;case 2D = L - C;PO(2) = PO(2) + Q*(L+2*C)*D^2/L^3;PO(3) = PO(3) + Q*C*D^2/L^2;PO(5) = PO(5) + Q-PO(2);PO(6) = PO(6) + Q*D*C^2/L^2;case 3D = L - C;PO(1) = PO(1) + Q*D/L;PO(4) = PO(4) + Q*C/L;case 4PO(2) = PO(2) + 7*Q*L/20;PO(3) = PO(3) + Q*L^2/20;PO(5) = PO(5) + 3*Q*L/20;PO(6) = PO(6) + Q*L^2/30;end子程序2function [KE, T] = ketb(dxy, E, A, IZ)DY = dxy(2,2) - dxy(1,2);DX = dxy(2,1) - dxy(1,1);L = sqrt(DX^2 + DY^2);S = DY/L;C = DX/L;a1 = IZ/L;a2 = a1/L;a3 = E/L;KE = a3*[A, 0, 0, -A, 0, 0 0, 12*a2, 6*a1, 0, -12*a2, 6*a10, 6*a1, 4*IZ, 0, -6*a1, 2*IZ-A, 0, 0, A, 0, 00, -12*a2, -6*a1, 0, 12*a2, -6*a10, 6*a1, 2*IZ, 0, -6*a1, 4*IZ]; t = [C, S, 0; -S, C, 0; 0, 0, 1];t1 = zeros(3,3);T = [t, t1; t1, t];end子程序function KZ = kztb(XY, ELB, EAIZ)[N,M] = size(XY);KZ =zeros(9, 9);[MB, m] = size(ELB);for k = 1 :2ii = ELB(k, 2);jj = ELB(k, 3);LTB = ELB(k, 4);dxy = [XY(ii,2), XY(ii,3)XY(jj,2), XY(jj,3)];E = EAIZ(LTB, 1);A = EAIZ(LTB, 2);IZ = EAIZ(LTB, 3);[KE, T] = ketb(dxy, E, A, IZ);CN = [3*ii-2, 3*ii-1, 3*ii, 3*jj-2, 3*jj-1, 3*jj];KE = (T')*KE*T;for i = 1:6for j = 1:6KZ(CN(i), CN(j)) = KZ(CN(i), CN(j)) + KE(i, j);endendend子程序function U = weiyi(KZ, P, SU)[LR, m] = size(SU);for k =1:LRi = SU(k, 1);P(i) = 1e10*SU(k, 2);KZ(i,i) = 1e10*KZ(i, i);endU = KZ\P;end子程序function P = ydx(XY, ELB, EAIZ, ELPQ)[N, m] = size(XY)P = zeros(3*N, 1)[LPQ, m] = size(ELPQ)for j = 1:LPQELPQ1 = ELPQ(j, :)k = ELPQ(j, 1)[PO, ii, jj] = dxjd(ELB, XY, ELPQ1)LTB = ELB(k, 4);dxy = [XY(ii,2), XY(ii,3)XY(jj,2), XY(jj, 3)];E = EAIZ(LTB, 1);A = EAIZ(LTB, 2);IZ = EAIZ(LTB, 3);[KE, T] = ketb(dxy, E, A, IZ );PO = (T')*PO;CN = [3*ii-2, 3*ii-1, 3*ii, 3*jj-2, 3*jj-1, 3*jj]for i = 1:6P(CN(i)) = P(CN(i)) + PO(i);endend后处理程序clc;[N, M] = size(XY);DU = zeros(N,4);for i = 1:NDU(i,1) = i;DU(i,2) = U(3*i-2);DU(i,3) = U(3*i-1);DU(i,4) = U(3*i);end[MB, M] = size(ELB);for i = 1:MBP1(i, 1) = i;end'节点号水平位移竖直位移转角'DU主程序>>GJINPUT;>>KZ = kztb(XY, ELB,EAIZ);>>P = ydx(XY, ELB, EAIZ, ELPQ);>>U = weiyi(KZ, P, SU);>>GJOUTPUT;运行结果：ans ='节点号水平位移竖直位移转角' DU =1.0000 -0.0000 -0.0000 -0.00002.0000 -0.0000 -0.0111 -0.01003.0000 -0.0231 -0.0111 -0.0125。
matlab桁架结构有限元计算
matlab桁架结构有限元计算摘要：一、引言二、MATLAB 在桁架结构有限元计算中的应用1.桁架结构的离散化与编号2.组装刚度矩阵3.求解器求解4.后处理三、MATLAB 桁架有限元分析实例1.空间桁架有限元分析2.基于MATLAB 的三维桁架有限元分析3.五杆桁架有限元分析四、用MATLAB 编写程序求解桁架结构内力问题1.静定结构绘制简图2.计算结构力学3.机动分析五、结论正文：一、引言在工程领域中，桁架结构由于其优越的力学性能和简单的结构形式，被广泛应用于桥梁、塔架等大型建筑结构中。
为了确保桁架结构的安全性和稳定性，对其进行有限元分析是非常必要的。
MATLAB 作为一种强大的数学软件，可以方便地进行桁架结构的有限元计算。
本文将介绍如何使用MATLAB 进行桁架结构有限元计算。
二、MATLAB 在桁架结构有限元计算中的应用1.桁架结构的离散化与编号在进行有限元分析之前，首先需要对桁架结构进行离散化处理，将连续的桁架杆件划分为有限个小杆件。
同时，为了方便后续计算，需要对各个杆件和节点进行编号。
2.组装刚度矩阵根据桁架结构的几何参数和材料性能，可以计算出各个杆件的刚度矩阵。
将这些刚度矩阵组装成一个总的刚度矩阵，用于描述整个桁架结构的刚度特性。
3.求解器求解利用MATLAB 的求解器，可以对桁架结构进行有限元求解。
求解器会根据刚度矩阵和施加的边界条件，计算出节点的位移和单元的应力。
4.后处理在求解完成后，需要对计算结果进行后处理。
这包括对计算结果的保存、可视化以及结果的验证等。
三、MATLAB 桁架有限元分析实例1.空间桁架有限元分析例如，可以针对一个空间桁架结构，使用MATLAB 进行有限元分析。
假设该桁架结构由铝制成的垂直和水平部分和钢制成的对角桁架构件组成。
结构承受荷载P，需要计算节点位移、单元应力以及支反力。
2.基于MATLAB 的三维桁架有限元分析还可以利用MATLAB 进行三维桁架有限元分析。
matlab有限元分析实例
matlab有限元分析实例问题描述：考虑一平面有界区域，设其边界为[。
我们求解泊松方程之狄利克雷边值问题。
问题的强形式为一椭圆型偏微分方程当之几何形状稍复杂时，一般无法求得其解析解。
我们可应用有限元法来求其数值解。
通常我们先写出该问题的抽象弱形式：求使得其中为检验函数为一适当索伯列夫空间（在本例中也是一希尔伯特空间）为一双线性型为一线性型。
其具体表达式为有限元空间离散我们采用最简单的二维单元离散单元即三节点线性三角形单元，其插值基函数（即形函数为一次多项式。
有限元之核心思想为：使用离散的函数空间来分片逼近连续的函数空间。
于是所求之近似解可写成基函数之线性组合其中为待求系数，常称为自由度。
将此近似解之表达式代入前述问题之弱形式，并取检验函数，可得写成矩阵形式，便成为常见之有限元方程由于历史原因，通常采用固体力学中的习惯命名：[公式] 为刚度矩阵，[公式] 为自由度向量，[公式] 为载荷向量。
由于空间已分片离散，上面的有限元方程只在各单元内部成立。
为了求解方便，通常我们将所有单元的有限元方程连立起来求解，于是需要将各单元之刚度矩阵组装成总体刚度矩阵。
求解器MATLAB 代码及简释本求解器之十行代码如下：function u=fem(nds,els,bcs)nnd=size(nds,1); u=zeros(nnd,1); K=zeros(nnd,nnd); f=zeros(nnd,1);for j=1:size(els,1)K(els(j,:),els(j,:))=K(els(j,:),els(j,:))+stima(nds(els(j,:),:));f(els(j,:))=f(els(j,:))+ones(3,1)*det([1,1,1;nds(els(j,:),:)'])/ 6;endfreends=setdiff(1:nnd,bcs);u(freends)=K(freends,freends)\f(freends);function stima=stima(vertices)ndim=size(vertices,2);J=[ones(1,ndim+1);vertices'];B=J\[zeros(1,ndim);eye(ndim)];stima=det(J).*B*B'/prod(1:ndim);输入数据格式：唯一的一个函数stima 用来计算各单元刚度矩阵。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
1.物理现象：这个对工程师来说是直观的物理现象和物理量，温
度多少度，载荷是多大等等。
通常来说，用户界面中呈现的、用户对工程问题进行设置时输入的都是此类信息。
2.数学方程：将物理现象翻译成相应的数学方程，例如流体对应
的是NS方程，传热对应的是传热方程等等；大部分描述这些现象的方程在空间上都是偏微分方程，偶尔也有ODE（如粒子轨迹、化学反应等）。
在这个层面，软件把物理现象“翻译”
为以解析式表示的数学模型。
3.数值模型：在定义了数学模型，并执行了网格剖分后，商业软
件会将数学模型离散化，利用有限元方法、边界元法、有限差分法、不连续伽辽金法等方法生成数值模型。
软件会组装并计算方程组雅可比矩阵，并利用求解器求解方程组。
这个层面的计算通常是隐藏在后台的，用户只能通过一些求解器的参数来干预求解。
有限元是一种数值求解偏微分方程的方法。
基本过程大致是设置形函数，离散，形成求解矩阵，数值解矩阵，后处理之类的。
MATLAB要把这些过程均自己实现，不过在数值求解矩阵时可以调用已有函数。
可以理解为MATLAB是一个通用的计算器，当然它的功能远不止如此。
而ANSYS之类的叫做通用有限元软件，针对不同行业已经将上述过程封装，前后处理也比较漂亮，甚至不太了解有限元理论的人也能算些简单的东西，当然结果可靠性又另说了。
比较两者，ANSYS之类的用起来容易得多，但灵活性不如MATLAB。
MATLAB用起来很困难，也有人做了一些模块，但大多数只能解决一些相对简单的问题。
对于大多数工程问题，以及某些领域的物理问题，一般都用通用有限元软件，这些软件还能添加一些函数块，用以解决一些需要额外设置的东西。
但是对于非常特殊的问题，以及一般性方程的有限元解，那只能用MATLAB或C，Fortran之类的了。