一、實驗題目
（1）用ritz-galerkin方法求解邊值問題
精確解（2）用有限元方法求解
二、實驗目的
運用matlab數學軟體編寫ritz-galerkin方法和有限元方法程式，進一步熟悉matlab的應用及掌握偏微分方程數值方法中ritz-galerkin方法和有限元方法，對各個方法求解精度進一步明確。
三、實驗原理
（1）ritz-galerkin方法
ritz方法
設給出二次泛函
（1）其中是hilbert空間v上的雙線性泛函，而且滿足對稱性、有界性、正定性；是v上的有界線性泛函。考慮一下變分問題：求滿足
2）設hilbert空間v是可分的，即v中有可數多個元素構成乙個稠密集。在v中取n個線性無關的元素，他們張成乙個n維子空間，即,或記成3）
上述元素稱為空間的基.
用代替v，在上求泛函的極小，即求滿足
4)顯然，原先的變分問題（2）與後面的變分問題（4）是不同的，他們的解也是不同的。如果v的子空間充分接近v，那麼，用代替v而得到的解也就充分接近u，從而把作為原變分問題的近似解，亦即把無窮維空間v中的極值問題近似為有限維空間中的極值問題，這就是ritz方法得基本思想。
問題（4）很容易化為乙個線性代數方程組問題。實際上，每乙個中的元素都可以用的線性組合表示，故，這時，由及的性質
，因此，是乙個以為變數的二次多項式。由於雙線性泛函是對稱、正定的，則二次項
是乙個以為變數的正定二次型，即矩陣是對稱正定的。因此，在點取得極小的充分必要條件是。
可以算出，上式即，即n維向量是線性方程組的解，其中。求出之後，即得近似變分問題(4)的解
總之，ritz法吧泛函的極小值問題化為多元二次函式的極小值問題，最後通過求解乙個線性代數方程組而得出變分問題的近似解。
galerkin方法
當雙線性泛函對稱時，變分問題（2）等價於變分方程：求使得（5）
galerkin方法是一種求解變分方程（5）的近似方法。假設v是可分的hilbert空間，仍然取有限維子空間為。考慮近似變分方程：求滿足6）
需要注意的是，近似變分方程（6）的解中，且檢驗函式v只要求對中的元素成立而不是對v中的一切元素成立。
這時，求解仍然化為求解乙個線性方程組。因為中的元素仍是通過的線性組合來表示，可以設，代入式（6）得
或者由的任意性，即得任意性，可得出線性代數方程組
由於雙線性泛函是對稱的，所以該線性方程組與ritz得到的線性方程組是相同的。
（1）有限元方法
考慮兩點邊值問題（p）
其中,且,， (記)，常數，為已知常數。
根據邊界條件，定義空間
構造空間v的有限維子空間線性無關的函式集合，使它成為空間的基。
引入空間上的雙線性泛函
和線性泛函
於是，有限元方法求關於問題（p）的解可以表示為
把區間剖分，然後作定義在區間、支集在和上的「山形」函式：
顯然線性無關，而且滿足，故由張成的子空間。
容易驗證
所以，如果給出節點上的未知函式值，那麼，在中的近似解就是。
我們的問題是：求使得
這相當於：求使得
定義矩陣
，， 則
可求出所以邊值問題的近似解為
四、實驗步驟
（1）ritz-galerkin方法
clear all
n=2;h=(1-0)/n;
syms x;
for i=1:n
fai(i)=x*(1-x)*(x^(i-1));
dfai(i)=diff(x*(1-x)*(x^(i-1)));
endfai;
for i=1:n
for j=1:n
fun=-dfai(i)*dfai(j)+fai(i)*fai(j);
a(i,j)=int(fun,x,0,1);
endfun=-x*fai(i);
f(i)=int(fun,x,0,1);
endc=inv(a)*f';
product=c.*fai';
c;h1=1/10;
sum=0;
for i=1:n
sum=sum+product(i);
endvu=;
for y=0:h1:1
v=subs(sum,x,y);
vu=[vu,v];
endy=0:h1:1;
yy=0:0.001:1;
u=sin(yy)/sin(1)-yy;
plot(yy,u,'-',y,vu,'o-');
u=vpa(u,5)
vu=vpa(vu,5)
iwanttoknow=abs(u-vu);
結果如下：
表1數值解與近似解對照表
圖1數值解與近似解對照圖
（2）有限元方法
clear all
syms x
n=5;h=1/(n+1);
for i=1:n
fai1(i)=1+(x-i*h)/h;
posai1(i)=diff(fai1(i));
fai2(i)=1+(i*h-x)/h;
posai2(i)=diff(fai2(i));
enda=zeros(n,n);
func111=posai1(1)^2+fai1(1)^2;func112=posai2(1)^2+fai2(1)^2;
a(1,1)=int(func111,x,0,h)+int(func112,x,h,2*h);
func12=posai2(1)*posai1(2)+fai2(1)*fai1(2);
a(1,2)=int(func12,x,h,2*h);
for i=2:n-1
func111=posai1(i)^2+fai1(i)^2;func112=posai2(i)^2+fai2(i)^2;
a(i,i)=int(func111,x,(i-1)*h,i*h)+int(func112,x,i*h,(i+1)*h);
func12=posai2(i)*posai1(i+1)+fai2(i)*fai1(i+1);
a(i,i+1)=int(func12,x,i*h,(i+1)*h);
func21=posai1(i)*posai2(i-1)+fai1(i)*fai2(i-1);
a(i,i-1)=int(func21,x,(i-1)*h,i*h);
endi=n;
func111=posai1(i)^2+fai1(i)^2;func112=posai2(i)^2+fai2(i)^2;
a(i,i)=int(func111,x,(i-1)*h,i*h)+int(func112,x,i*h,(i+1)*h);
func21=posai1(i)*posai2(i-1)+fai1(i)*fai2(i-1);
a(i,i-1)=int(func21,x,(i-1)*h,i*h);
for i=1:n
funcb1=fai1(i)*(1+pi^2)*sin(pi*x); funcb2=fai2(i)*(1+pi^2)*sin(pi*x);
b(i)=int(funcb1,x,(i-1)*h,i*h)+int(funcb2,x,i*h,(i+1)*h);
endb=b';
a;u=inv(a)*b;
u=vpa(u,5);
u=u';
x=h:h:1-h;
xx=0:0.001:1;
y=sin(pi*xx);
plot(xx,y,'-',x,u,'o-')
結果如下：
表2數值解與近似解對照表
圖2數值解與近似解對照圖
微分方程數值解法實驗報告
for k 2 n yp y k 1,h feval f,x k 1 y k 1,yc y k 1,h feval f,x k yp y k,yp yc 2 end 3 runge kutta method function tout,yout runge kutta f,ab,y0,h tout ...
常微分方程數值解實驗報告
學院 數學與資訊科學 專業 資訊與計算科學 姓名 鄭思義 學號 201216524 課程 常微分方程數值解 實驗一 常微分方程的數值解法 1 分別用euler法 改進的euler法 預報校正格式 和s k法求解初值問題。h 0.1 並與真解作比較。1.1實驗 尤拉法 function x,y nae...
常微分方程數值解實驗報告
2012年10月17日 實驗目的 1 通過用matlab程式設計運用euler方法及其改進方法求解常微分方程初值問題，更進一步掌握常微分方程及其數值解法課程的理論內容，加深對數值解法的理解。2 熟悉matlab程式設計環境。實驗內容 1.實驗題目 運用euler方法及其改進方法求解常微分方程初值問題...