有限差分法c++實現
幾種eluer演算法
在matlab 裡面寫乙個程式:要求用 隱式尤拉法(backward euler) 去解決常微分方程。
下面是兩個例題,
1. x' = - 2x ; 還給出準確值是: x= e^(-2t) 要求: 求出t=0, 這是乙個初始值,然後算出在區間[0,5] 的值。
還給出 步長(step size) h=0.1
2. x' = xsint - 2x^2; 給出初始條件 x(0) = 0
要求: 求出在區間[0,1] 的值,比較 h=0.01 和 h=0.001 在 t = 1 這裡。
1.新建乙個m檔案,編寫隱式euler法的程式:
function [x,y]=implicit_euler(odefun,xspan,y0,h,varargin)
% 隱式euler公式求解常微分方程
% 輸入引數:
% ---odefun:微分方程的函式描述
% ---xspan:求解區間[x0,xn]
% ---y0:初始條件
% ---h:迭代步長
% ---p1,p2,…:odefun函式的附加引數
% 輸出引數:
% ---x:返回的節點,即x=xspan(1):h:xspan(2)
% ---y:微分方程的數值解
x=xspan(1):h:xspan(2);
y(1)=y0;
for k=1:length(x)-1
z0=y(k)+h*feval(odefun,x(k),y(k),varargin);
z1=inf;
while abs(z1-z0)>1e-4
z1=y(k)+h*feval(odefun,x(k+1),z0,varargin);
z0=z1;
end
y(k+1)=z1;
end
x=x;y=y;
2.在命令視窗直接呼叫上面的程式,求解常微分問題:
(1)
f=@(t,x)-2*x;
[t,x]=implicit_euler(f,[0 5],1,0.1);
t1=0:0.1:5;
x1=exp(-2*t);
plot(t1,x1)
hold on
plot(t,x,'k.-','markersize',16)
legend('解析解','隱式euler求解結果')
xlabel('t');ylabel('x');
(2)此題給出的初值好像有問題吧,x0=0的話,求解的結果都是為0,所以我改用x0=1求解試了一下:
f=@(t,x)x*sin(t)-2*x^2;
[t,x]=implicit_euler(f,[0 1],1,0.01);
[t1,x1]=implicit_euler(f,[0 1],1,0.001);
e=x1(end)-x(end)
>>e =
-0.0029
plot(t,x,'r:')
hold on
plot(t1,x1,'g--')
xlabel('t');ylabel('x')
legend('積分步長為0.01','積分步長為0.001')
利用matlab實現用:
向前euler方法
向後euler方法
改進的euler方法
梯形公式
求方程 y'=y-2x/y, 0
有限差分法MATLAB程式
設有乙個長直接地金屬矩形槽,長a 40,寬b 20,其側壁與底面電位均為零,頂蓋電位為100v 相對值 求槽內電位分布。利用高斯迭代求解 如下 相鄰兩次迭代值最大允許誤差為0.001 a zeros 21,41 a 1,100 b zeros 19,39 c eye 19,39 count 1 d ...
有限元法 有限差分法 有限體積法
有限元法也叫有限單元法 finite element method,fem 是隨著電子計算機的發展而迅速發展起來的一種彈性力學問題的數值求解方法。五十年代初,它首先應用於連續體力學領域 飛機結構靜 動態特性分析中,用以求得結構的變形 應力 固有頻率以及振型。由於這種方法的有效性,有限單元法的應用已從...
有限元法,有限差分法和有限體積法的區別
原文 有限差分方法 fdm 是計算機數值模擬最早採用的方法,至今仍被廣泛運用。該方法將 求解域劃分為差分網格,用有限個網格節點代替連續的求解域。有限差分法以taylor級 數展開等方法,把控制方程中的導數用網格節點上的函式值的差商代替進行離散,從而 建立以網格節點上的值為未知數的代數方程組。該方法是...