Symbolic Computing and Innovation of Engineering Education / Researching
符號計算 与 工程教育/研究創新

Symbolic Computing and Innovation of Engineering Education / Researching
符號計算 与 工程教育/研究創新

Symbolic Computing - the IDE of the next generation
符號計算 - 下一代的 集成開發環境

Symbolic Computing - the IDE of the next generation
符號計算 - 下一代的 集成開發環境

What does the IDE of the next generation look like? Is it similar to C++ and Fortran?
下一代的集成開發環境是怎樣的？和我們熟悉的C++、Fortran、Java有哪些顯著的區別？

Now you are just watching the next generation IDE, which is purely symbolic, and compatible with C++ and Fortran.
您看到的正是下一代的集成開發環境。它完全基於符號系統，同時兼容 C++ 和 Fortran。

Grid[{{CForm[Expand[(a+b)^3]]}}]

Grid[{{FortranForm[Expand[(a+b)^3]]}}]

Which C++ and Fortran can not do is:
C++ 和 Fortran 幹不了的活：

Grid[{{Expand[(a+b)^3]}}]

and more:
還有：

Style[Manipulate[Grid[{{Expand[(a+b)^n]}}],{n,1,100,1}],FontColor->Black]

Symbolic Computing - the IDE of the next generation
符號計算 - 下一代的 集成開發環境

What does the next generation IDE look like? Is it similar to C++ and Fortran?
下一代的集成開發環境是怎樣的？和我們熟悉的C++、Fortran、Java有哪些顯著的區別？

and, the most important thing in engineering computing - The precision!
還有，在工程計算中最重要的東西 － 精度！

Style[Manipulate[Grid[{{"\nPi in 64 bit Symbolic Computing: \n","\nPi in 64 bit Numerical Computing: \n"},{N[Pi,n],N[Pi,308]}}],{n,308,10000,1}],FontColor->Black]

Symbolic Computing - the IDE of the next generation
符號計算 - 下一代的 集成開發環境

How can we change our researching and teaching habits to fit the high-efficiency researching workflow in the future?
我們應該如何改變我們習慣的研究及教學方式，以適應未來高效的研究流程？

1. Use symbolic computing TODAY
   用符號計算
    
2. Pay close attention to the original mathematical principles
   關注問題背後數學定理
    
3. Make the new workflow
   制定新流程

Engineering Education/Researching is Transforming
工程教育/研究正在轉型

The Traditional Workflow of Engineering Education:
傳統的工程教學流程:

1. The teacher performs the mathematical principles of engineering problems
   老師在課堂上展示工程算法的基本數學原理
   
2. The teacher use examples to show how to solve simplified engineering problems with mathematical principles
   老師利用例題闡述如何利用基本原理解決簡化了的工程問題

3. Since the real engineering problems are so complicated, they are impossible to be solved by hand, the teacher will told the students to use commercial softwares to solve real engineering problems
   由於現實中的工程問題十分複雜，利用基本數學原理手工推導求解幾乎不可能實現，老師一般會教學生如何使用成熟的商業軟件解決現實中的工程問題

Consequences:
後果:

1. Students have no experiences of solving engineering problems by using original mathematical principles. They do not fully understand the principles behind the engineering problems.
   學生沒有利用基本數學原理解決現實工程問題的經驗，在基本數學原理与現實工程問題之間缺乏貫通。學生對現實工程問題背後的數學原理理解不夠深刻。

2. Since the students have been told to use commercial softwares to solve engineering problems, they will just remember the workflows of commercial softwares, but forget the mathematical principles.
   由於學生被告知利用成熟的黑箱式商業軟件解決現實中的工程問題，久而久之學生只能記住商業軟件的操作流程，而忘記其背後深刻的數學原理。

3. Since the students can not fully understand the mathematical principles behind engineering problems, they can not innovate at all.
   由於學生不能充分理解工程問題背後的數學原理，學生便不可能在前人的基礎上有所創新。

Engineering Education/Researching is Transforming
工程教育/研究正在轉型

The New Workflow of Engineering Education:
新的工程教學流程:

1. The teacher performs the mathematical principles of engineering problems
   老師在課堂上展示工程算法的基本數學原理
   
2. The teacher use examples to show how to solve simplified engineering problems with mathematical principles
   老師利用例題闡述如何利用基本原理解決簡化了的工程問題

5. The teacher teach the students how to use commercial softwares to solve real engineering problems, the students will understand the mathematical principles behind these commercial softwares
   老師教會學生如何使用成熟的商業軟件解決現實中的工程問題，同時讓學生理解商業軟件背後的數學原理

Result:
結果:

Sample - How to use Symbolic IDE to innovate with Weighted Residual Method & Finite Element Method
範例 － 如何利用符號計算創新加權殘值法与有限單元法

The Application of Weighted Residual Method ( WRM ) & Finite Element Method ( FEM )
加權殘值法与有限單元法的應用

Engineering Cases:
工程案例:

Analysis of the vibration of the solar panel on a spaceship:
太空船太陽能電池板的共振分析

FEM Model:
有限元模型:

The vibration of the solar panel - 64Hz
太空船太陽能電池板的共振情況 - 64Hz

Import["~/Documents/Solar.avi","Animation"]

WRM and FEM have been widely used in aerospace, automobile, ship building and other manufacturing areas. Computer-Aided-Engineering (CAE), which is based on these two methods, has made great contributions to the productivity of mankind.
加權殘值法 及 有限單元法 已被大規模應用於 航空航天、汽車、船舶等製造領域。以此兩種方法為基礎發展而來的 計算機輔助工程（CAE） 行業為人類生產力的發展做出了卓越的貢獻。

But, CAE is just the beginning.....
但是，CAE 僅僅是個開始.......

The Application of Weighted Residual Method ( WRM ) & Finite Element Method ( FEM )
加權殘值法与有限單元法的應用

The modern workflow of Product Development:	
現代產品研發的工作流程:

Conceptual Design -> CAD -> CAE -> Physical Prototype
概念設計 -> 計算機輔助設計 -> 計算機輔助工程 -> 物理原形

Data Driven vs Equation Driven
數據驅動 vs 方程驅動

To solve engineering problems, we just have to input some parameters into the commercial software. The mathematical equations behind the software are fixed. For example, the original equation of MSC.Nastran 101 Solver is  (([δ])^T)∫(([B])^T)[D][B]dA [δ] = (([δ])^T)[F]. This kind of method is called "Data Driven". Most CAE softwares are Data Driven Softwares.

With symbolic IDE, we can break the limits of "Data Driven". We can use our mathematical methods - WRM and FEM - onto every kind of equations. This is called "Equation Driven".

商業軟件只需我們給定參數便可求解工程問題，但其背後的數學方程的形式是不變的。例如 MSC.Nastran 101求解器的主要數學方程便是 (([δ])^T)∫(([B])^T)[D][B]dA [δ] = (([δ])^T)[F]。這種只改參數不改方程的形式稱為“數據驅動”。大多數 CAE 軟件是數據驅動軟件

符號開發環境讓我們有可能打破“數據驅動”的局限性，將數學方法作用到任意的方程或者函數上，這種形式稱為“函數驅動”。

Advantages of Equation Driven:

Not only the engineering problems that can be solved in Data-Driven Softwares, but also the problems which can not be covered by Data-Driven Softwares, can be solved by Equation-Driven Softwares

方程驅動的好處：

方程驅動的軟件不僅能夠解決數據驅動軟件能夠解決的工程問題，更能夠解決當前數據驅動軟件不能夠涵蓋的問題。

To develop Equation-Driven Softwares, we need:

1. The knowledge of general mathematical principles

2. Pure Symbolic Development Environment

要開發方程驅動的軟件，我們需要：

1. 通用的數學知識

2. 純符號的開發環境

Weighted Residuals Method - 加權殘值法

First, we have to define the differential equation in Mathematica.
首先，我們需要在 Mathematica 中定義需要求解的微分方程。

Sample in Mathematica - How to define a Function :
Mathematica 範例 - 如何定義一個函數：

DifferentialFunction[deflection_]:=(∂_x)(∂_x)deflection;
Grid[{{"The Equation is: ",DifferentialFunction[v[x]],"==",(ℽ/ElaModulus)}}]

微分算子[撓度_]:=(∂_x)(∂_x)撓度;
Grid[{{"方程是: ",微分算子[v[x]],"==",(ℽ/(彈性模量))}}]

Mathematica can solve differential equations directly, by using DSolve. But that's not our point today.
雖然 Mathematica 可以使用 DSolve 直接求解微分方程，但這不是今天的重點。

Style[DSolve[DifferentialFunction[v[x]]==(ℽ/ElaModulus),v[x],x]//MatrixForm]

Weighted Residuals Method - 加權殘值法

What is the Residual?
殘值是甚麼？

DifferentialFunction[deflection_]:=(∂_x)(∂_x)deflection;
Grid[{{"The Residual is: ",DifferentialFunction[(OverTilde[v])[x]] - (ℽ/ElaModulus),"≠", 0}}]

Weighted Residuals Method - 加權殘值法

(∫_A)Residual*Weight[i]dA=0,     i = 1, 2, 3, ⋯⋯, n

Grid[{{∫((OverTilde[v])''[x]-(ℽ/ElaModulus))W[i]\[DifferentialD]x,"==",0}}]

Weighted Residuals Method - 加權殘值法

Collocation Method - 配點法

Weight[i](x)=Dirac(x-(Subscript[x, i])), i = 1, 2, 3, ⋯⋯, n

Grid[{{(OverTilde[v])''[(Subscript[x, i])]-(ℽ/ElaModulus),"==",0}}]

Weighted Residuals Method - 加權殘值法

Sub-domain Method - 子域法

(∫_A)Residual*Weight[i]dA=∑(∫_(Subscript[A, i]))Residuald(Subscript[A, i]), i = 1, 2, 3, ⋯⋯, n

Grid[{{"(∑_(i = 1)^n)∫(((OverTilde[v])^(′′))[(Subscript[x, i])]-(ℽ/ElaModulus))\[DifferentialD](Subscript[x, i])","==",0}}]

Weighted Residuals Method - 加權殘值法

Least Squares Method - 最小二乘法

Weight[i]=(∂_(Subscript[α, i]))Residual, i = 1, 2, 3, ⋯⋯, n

Style["∫(((OverTilde[v])^(′′))[x]-(ℽ/ElaModulus))(∂_(Subscript[α, i]))((OverTilde[v])''[x]-(ℽ/ElaModulus))\[DifferentialD](Subscript[x, i]) = 0"]

Weighted Residuals Method - 加權殘值法

Galerkin Method - 迦遼金法

Weight[i]=(∂_(Subscript[α, i]))(OverTilde[ U]), i = 1, 2, 3, ⋯⋯, n

Style["∫(((OverTilde[v])^(′′))[x]-(ℽ/ElaModulus))(∂_(Subscript[α, i]))((OverTilde[v])''[x])\[DifferentialD](Subscript[x, i]) = 0"]

Weighted Residuals Method - 加權殘值法

Example - 例子：

Use all four methods above to solove the differential equation below - 利用上述四種方法求解下列微分方程：

(∂_x)(∂_x)( U (x) )+U =1; U (0)=1; U (1)=0.

Basis Function: (x^i) ,  i=0, 1, 2,⋯⋯, n 
基函數為：(x^i) ,  i=0, 1, 2,⋯⋯, n

DifferentialFunction[var_]:=(∂_x)(∂_x)var+var;
Residual[U_]:=DifferentialFunction[U]-1;
U[x_]:=(Subscript[α, 0])+(Subscript[α, 1])x+(Subscript[α, 2])(x^2)+(Subscript[α, 3])(x^3)+(Subscript[α, 4])(x^4);
Grid[{{"The Residual is: ",Residual[Evaluate[U[x]]]}}]

According to the boundary conditions, we have two equations below: - 由邊界條件，我們得到以下兩個方程：

Grid[{{U[0]," = ",1},{U[1]," = ",0}}]

Three more equations are needed. - 還需要三個方程式
We will use different methods to choose the weights. - 我們將使用不同的方法取權

ResidualC[x_]:=Evaluate[Residual[U[x]]];
Grid[{{ResidualC[0.25]," = ",0},{ResidualC[0.5]," = ",0},{ResidualC[0.75]," = ",0}}]
(Subscript[Result, 1])=NSolve[{ResidualC[0.25]==0,ResidualC[0.5]==0,ResidualC[0.75]==0,(Subscript[α, 0])== 1,(Subscript[α, 0])+(Subscript[α, 1])+(Subscript[α, 2])+(Subscript[α, 3])+(Subscript[α, 4])==0},{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}];
Grid[{{Transpose[{{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}}]//MatrixForm," = ",Transpose[{Table[(Subscript[Result, 1])[[1,i,2]],{i,5}]}]//MatrixForm}}]
collocation[x_]:=Table[(Subscript[Result, 1])[[1,i,2]],{i,5}].{1,x,(x^2),(x^3),(x^4)};

Grid[{{(∫_0^0.33)Residual[U[x]]\[DifferentialD]x," = ",0},{(∫_0.33^0.66)Residual[U[x]]\[DifferentialD]x," = ",0},{(∫_0.66^1)Residual[U[x]]\[DifferentialD]x," = ",0}}]
(Subscript[Result, 2])=NSolve[{(∫_0^0.33)Residual[U[x]]\[DifferentialD]x==0,(∫_0.33^0.66)Residual[U[x]]\[DifferentialD]x==0,(∫_0.66^1)Residual[U[x]]\[DifferentialD]x==0,(Subscript[α, 0])== 1,(Subscript[α, 0])+(Subscript[α, 1])+(Subscript[α, 2])+(Subscript[α, 3])+(Subscript[α, 4])==0},{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}];
Grid[{{Transpose[{{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}}]//MatrixForm," = ",Transpose[{Table[(Subscript[Result, 2])[[1,i,2]],{i,5}]}]//MatrixForm}}]
subdomain[x_]:=Table[(Subscript[Result, 2])[[1,i,2]],{i,5}].{1,x,(x^2),(x^3),(x^4)};

(Subscript[weight, 1])=(∂_(Subscript[α, 2]))Residual[U[x]];(Subscript[weight, 2])=(∂_(Subscript[α, 3]))Residual[U[x]];(Subscript[weight, 3])=(∂_(Subscript[α, 4]))Residual[U[x]];
Grid[{{(∫_0^1)Residual[U[x]](Subscript[weight, 1])\[DifferentialD]x," = ",0},{(∫_0^1)Residual[U[x]](Subscript[weight, 2])\[DifferentialD]x," = ",0},{(∫_0^1)Residual[U[x]](Subscript[weight, 3])\[DifferentialD]x," = ",0}}]
(Subscript[Result, 3])=NSolve[{(∫_0^1)Residual[U[x]](Subscript[weight, 1])\[DifferentialD]x==0,(∫_0^1)Residual[U[x]](Subscript[weight, 2])\[DifferentialD]x==0,(∫_0^1)Residual[U[x]](Subscript[weight, 3])\[DifferentialD]x==0,(Subscript[α, 0])== 1,(Subscript[α, 0])+(Subscript[α, 1])+(Subscript[α, 2])+(Subscript[α, 3])+(Subscript[α, 4])==0},{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}];
Grid[{{Transpose[{{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}}]//MatrixForm," = ",Transpose[{Table[(Subscript[Result, 3])[[1,i,2]],{i,5}]}]//MatrixForm}}]
leastsquares[x_]:=Table[(Subscript[Result, 3])[[1,i,2]],{i,5}].{1,x,(x^2),(x^3),(x^4)};

(Subscript[weight, 1])=(∂_(Subscript[α, 2]))U[x];(Subscript[weight, 2])=(∂_(Subscript[α, 3]))U[x];(Subscript[weight, 3])=(∂_(Subscript[α, 4]))U[x];
Grid[{{(∫_0^1)Residual[U[x]](Subscript[weight, 1])\[DifferentialD]x," = ",0},{(∫_0^1)Residual[U[x]](Subscript[weight, 2])\[DifferentialD]x," = ",0},{(∫_0^1)Residual[U[x]](Subscript[weight, 3])\[DifferentialD]x," = ",0}}]
(Subscript[Result, 4])=NSolve[{(∫_0^1)Residual[U[x]](Subscript[weight, 1])\[DifferentialD]x==0,(∫_0^1)Residual[U[x]](Subscript[weight, 2])\[DifferentialD]x==0,(∫_0^1)Residual[U[x]](Subscript[weight, 3])\[DifferentialD]x==0,(Subscript[α, 0])== 1,(Subscript[α, 0])+(Subscript[α, 1])+(Subscript[α, 2])+(Subscript[α, 3])+(Subscript[α, 4])==0},{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}];
Grid[{{Transpose[{{(Subscript[α, 0]),(Subscript[α, 1]),(Subscript[α, 2]),(Subscript[α, 3]),(Subscript[α, 4])}}]//MatrixForm," = ",Transpose[{Table[(Subscript[Result, 4])[[1,i,2]],{i,5}]}]//MatrixForm}}]
galerkin[x_]:=(Table[(Subscript[Result, 4])[[1,i,2]],{i,5}].{1,x,(x^2),(x^3),(x^4)});

ps={Red,Blue,Thick,Dashed,Directive[Red,Thick]};
Plot[{collocation[x],subdomain[x],leastsquares[x],galerkin[x]},{x,0,1},PlotStyle->ps]
Plot[{collocation[x],subdomain[x],leastsquares[x],galerkin[x]},{x,0.49,0.51},PlotStyle->ps]

Weighted Residuals Method - 加權殘值法

Interactive Time - 互動時間

What about the basis function? Can we choose any basis function without limitation?
形函數呢？我們能否隨意選取形函數？

In Weighted Residuals Method, different basis functions can cause widely divergent results.
在加權殘值法中，不同的形函數導致的結果將大相徑庭。
It's very convenience to check the results caused by different basis functions in Mathematica interactively just by using "Manipulate".
在 Mathematica 中，我們能夠非常方便的以互動形式檢驗不同形函數所導致的結果 － 只需調用"Manipulate"。

Example - 例子

Differential Equation - 微分方程：(∂_x)(∂_x)v(x) = 1

Differential2(*定義待解微分方程的名稱為“Differential2”*)[var_]:=(∂_x)(∂_x)var;

Basis Functions - 形函數：

(Subscript[Sin, 1])[var_]:=Sin[var];
XBased[var_]:=var(var-1);
XTerm[var_,n_]:=var^(n);

Collocation Method - 配點法：

WRM[shapefunction_,function_,force_,pointN_,i_]:=
(
shapeT[i]=Table[shapefunction,{N,pointN}];
residualT[i]=function[shapeT[i]];
matrixK[i]=Table[residualT[i],{x,1/(pointN+1),1-1/(pointN+1),1/(pointN+1)}];
rightT[i]=Table[1,{N,pointN}];
parameter[i]=Table[(Subscript[a, j]),{j,pointN}];
showBase[i]=shapeT[i];
symbolicalresult[i]=shapeT[i].parameter[i];
numericalparameter[i]=N[Inverse[matrixK[i]].rightT[i]];
resultT[i]=shapeT[i].numericalparameter[i];
);

Maniputale - 互動操作

Style[Manipulate[
Grid[{
{Grid[{{"\nCase 1: \n Shape Function: "},{basefunction}}]},
{
WRM[basefunction,Differential2,force,pointN,1];Grid[{{"Shape Functions:"},{showBase[1]//TableForm},{"\nParameters:"},{Grid[{{parameter[1]//TableForm,Grid[Transpose[{Table["=",{i,pointN}]}]],numericalparameter[1]//TableForm}}]}}],
Plot[{resultT[1],0.5x^2-0.5x},{x,0,1},Filling->Axis,PlotRange->{{0,1},{-0.2,0.2}},ImageSize->400]
},
{"",
Grid[{{"Equation Set:",N[matrixK[1]]//MatrixForm,parameter[1]//MatrixForm,Grid[Transpose[{Table["=",{i,1}]}]],rightT[1]//MatrixForm}}]
}
}],
{{pointN,2,"Number of Terms:"},Range[2,4,2],ControlType->Setter},
{{basefunction,XBased[x] Sin[N π x],"Shape Function:"},{XBased[x] Sin[N π x],Sin[N π x],Sin[N π x]XTerm[x,N],Cos[N π x]},ControlType->PopupMenu}
],FontColor->Black]

Weighted Residuals Method - 加權殘值法

Let's reverse our mind - 把事情反過來想

What we have done on the computer is the very task that the computer have to do for us!
我們在計算機上實現的東西正是計算機需要幫我們完成的任務！

In Weighted Residuals Method, we need the computer to pick the best basis function out for us.
在加權殘值法中，我們需要計算機為我們決定哪一個形函數最好。

Feel the power of symbolic computing!
感受符號計算的力量！

Example - 例子

Which basis function is better?
哪一個形函數更好？

Define the basis function (shape function) library - 定義形函數庫:

shapefunctionc={XBased[x]Sin[N π x],Sin[N π x],Cos[N π x]};

Ask the computer to find the best basis function out from the basis function library - 讓計算機在形函數庫中找出最好的形函數:

forcec=1;
pointNc=3;
wrmc[i_]:=(
shapeTc=Table[shapefunctionc[[i]],{N,pointNc}];
residualTc=Differential2[shapeTc];
matrixKc=Table[residualTc,{x,1/(pointNc+1),1-1/(pointNc+1),1/(pointNc+1)}];
rightTc=Table[forcec,{N,pointNc}];
varc=N[Inverse[matrixKc].rightTc];
resultTc=shapeTc.varc;
resultc[x_]:=Evaluate[resultTc];
resultc[0.5]
)
errorc=10;
For[i=1,i≤2,i++,If[Abs[wrmc[i]-(-0.12)]<errorc,ii=i;errorc=Abs[wrmc[i]-(-0.12)]]];
Grid[{{"Shape Function","'",shapefunctionc[[ii]],"'","is better!"},{"形函數","'",shapefunctionc[[ii]],"'","更好!"}}]

Example - 例子
Integration by part - 分步積分

∫(a (x^4)+b x)×(((∂^3)(c x+d))/(∂(x^3)))dx - zero or a single constant is no use for engineering computing. 零或者單一一個常數對於工程計算來說是沒有用的。

∫(a (x^4)+b x)D[c x+d,{x,3}]\[DifferentialD]x

We need to see the whole process of integration.
我們需要看到積分整個過程。

AbstractIntegrateD[function_,{x_,order_}]:=D[function,{x,order-1}];(*∫((\[DifferentialD] function)/(\[DifferentialD]x))\[DifferentialD]x = function*)
IntegratebyPart[partA_,partB_,{x_,order_}]:=partA.AbstractIntegrateD[partB,{x,order}]-Integrate[D[partA,{x,1}].AbstractIntegrateD[partB,{x,order}],x];
SuperIntegratebyPart[partA_,partB_,{x_,order_}]:=(n=0;
partI=partA;
partII=partB;
orderN=order;
While[AbstractIntegrateD[partII,{x,orderN}]===0&&partI=!=0,
n++;
partI=D[partI,x];
orderN--;
];
(-1)^n*IntegratebyPart[partI,partII,{x,orderN}]
);
Grid[{{"Result: ",SuperIntegratebyPart[a x^4+b x,c x+d,{x,3}]}}]

Why is zero or a single constant no use for us?
為甚麼 0 或者單一的常數沒有用？

The answer is: the element stiffness matrix, which is the core element of Finite Element Method, is hidden in the process of Integration-by-part.
答案是：有限單元法的核心－單元剛度矩陣，就隱藏在分步積分的過程中。

Finite Element Method - 有限單元法

Weighted Residual Finite Element Method - 加權殘值有限元法

Weight[i]=(∂_x)N [x]

Example - 例子

In this example, we'll show that we can use the computer to automatically derive the element stiffness matrix, just by giving the basis functions to the computer. ( This kind of task can hardly be done in ANSYS or MSC.Patran )
本例將演示只要給出形函數，我們便可以利用計算機自動推導單元剛度矩陣。( 此類任務在 ANSYS 或者 MSC.Patran 中難以完成 )

Differential Equation to be solved - 待解微分方程：   (∂_x)(∂_x)( U ) = f ( x )
Elements to be used - 使用單元類型：BAR - 杆單元 (two nodes in each element, one dof on each node - 每單元兩個節點，每節點一個自由度)

(OverTilde[ U]) = (Subscript[α, 0])+ (Subscript[α, 1])x;   (OverTilde[ U]) = (Subscript[N, 1])[x](Subscript[U, p])+ (Subscript[N, 2])[x](Subscript[U, p])

Coor={{1,(Subscript[X, p])},{1,(Subscript[X, q])}};
CoorBasisFunction={{1,x}};
Deflection={{(Subscript[U, p])},{(Subscript[U, q])}};
BasisFunction=FullSimplify[CoorBasisFunction.Inverse[Coor]];
Grid[{{"(OverTilde[U]) =",BasisFunction//MatrixForm,Deflection//MatrixForm}}]

Residual (殘值) =  (∂_x)(∂_x)( (OverTilde[ U]) ) - f ( x )；Weight(權值)=(∂_x)N [x]

(∫_A)Residual*WeightdA=0,

ResidualFEM=FullSimplify[D[BasisFunction.Deflection,{x,2}]-f];
Grid[{{"Residual: ",ResidualFEM//MatrixForm}}]

In this case (∂_x)(∂_x)(OverTilde[ U]) = 0, So we must use Integration by Part to conquer this problem.
在此例中，(∂_x)(∂_x)(OverTilde[ U])=0，因此必須使用分步積分來解決以上問題。

BasisFunctionTrans=Transpose[BasisFunction];
SuperIntegratebyPart[BasisFunctionTrans,BasisFunction,{x,2}];
ElementStiffness=FullSimplify[Integrate[D[partI,{x,1}].AbstractIntegrateD[partII,{x,orderN}],{x,(Subscript[X, p]),(Subscript[X, q])}]];
Grid[{{"element stiffness matrix is: ",ElementStiffness//MatrixForm}}]

Finite Element Equation should be:
有限元方程為：

Grid[{{ElementStiffness//MatrixForm," ",Deflection//MatrixForm," = ",{{ForceI},{ForceII}}//MatrixForm}}]

Finite Element Method - 有限單元法

What about meshing? - 划網格呢？

Plot3D[Sin[x^2+y^2],{x,-2,2},{y,-3,3},Mesh->All]

Thanks to Mathematica's powerful "Plot" function, it can mesh any surface into triangle elements. But, there's a little problem - we can not get the mesh data directly.
感謝Mathematica中強大的"Plot"命令，它能夠將任何面體劃分為三角形單元。但是，有一個小問題 － 我們無法直接獲得網格的數據。

To solve this problem, we must think about one basic concepet in computer science - DATA Strcture.
要解決這個問題，我們必須考慮使用計算機科學中的一個基本概念 － 數據結構。

plotsimple=Plot3D[Sin[x^2+y^2],{x,-2,2},{y,-3,3},Mesh->All];
datasimple=InputForm[plotsimple];

By using "InputForm", we can see the source code of every function.
通過使用"InputForm"，我們可以看到任何函數的源代碼。

"TreeForm" can show the tree structure of the source code.
"TreeForm"可以將源代碼的樹狀結構顯示出來。

So, we can easily draw out the meshing data.
因此，我們可以方便地將網格信息抽出來。

TreeForm[datasimple,VertexLabeling->Automatic,ImageSize->350,AspectRatio->1,PlotStyle->Directive[PointSize[Large],Red]]

By using "Manipulate", we can even design our own meshing GUI.
通過使用"Manipulate"，我們甚至可以設計自己的網格劃分界面。

Example - 例子

Style[Manipulate[
Grid[
{
{
plot=Plot3D[Sin[x^c+y^d],{x,-3,3},{y,-2,2},Mesh->None,ImageSize->300,AspectRatio->1,ViewPoint->Top,Axes->None],
data=InputForm[plot];
Grid[{{
"Nodes:","","Element Number:"
},
{
nodeN=Dynamic[Length[data[[1,1,1]]]],"",elemN=Dynamic[Length[data[[1,1,2,1,1,2,1,1,1]]]]
}}]
},
{
Show[Graphics3D[{EdgeForm[{Thick,Red}],GraphicsComplex[Dynamic[data[[1,1]][[1]]],Polygon[Dynamic[data[[1,1]][[2,1]][[1,2]][[1]][[1]][[1]][[elementnumber]]]]]}],Graphics3D[GraphicsComplex[Dynamic[data[[1,1]][[1]]],Polygon[Dynamic[data[[1,1]][[2,1]][[1,2]][[1]][[1]][[1]]]]]
],ImageSize->300,AspectRatio->1,ViewPoint->Top,Axes->None],
Grid[{
{
{{"Coordinates of Element No.",Dynamic[elementnumber]}}//TableForm
},
{
Grid[{{coor1=Dynamic[data[[1,1,1]][[data[[1,1,2,1,1,2,1,1,1,elementnumber,1]]]]]},
{coor2=Dynamic[data[[1,1,1]][[data[[1,1,2,1,1,2,1,1,1,elementnumber,2]]]]]},
{coor3=Dynamic[data[[1,1,1]][[data[[1,1,2,1,1,2,1,1,1,elementnumber,3]]]]]}
}]
}
}]
}
}
,Alignment->{Center}]
,{{c,0,"X Order"},0,10,1,ContinuousAction->False},{{d,0,"Y Order"},0,10,1,ContinuousAction->False},{{elementnumber,1,"Element Number"},1,elemN,1},ControlPlacement->Bottom],FontColor->Black]

Finite Element Method - 有限單元法

The researching workflow can be largely impoved, if we can master the symbolic computing.
掌握符號計算，研究流程將會大幅改進。

The original researching workflow is:
原先的研究流程是：



The new researching workflow will be:
新的研究流程將是：

The Application of Weighted Residual Method ( WRM ) & Finite Element Method ( FEM )
加權殘值法与有限單元法的應用

Mathematics - based engineering knowledge sharing:
基於數學的工程知識共享體系:

Symbolic Computing and Innovation of Engineering Education / Researching
符號計算 与 工程教育/研究創新

Symbolic Computing and Innovation of Engineering Education / Researching
符號計算 与 工程教育/研究創新