偏微分方程数值算法和matlab实验报告.docx
- 文档编号:15141320
- 上传时间:2023-07-01
- 格式:DOCX
- 页数:9
- 大小:71.35KB
偏微分方程数值算法和matlab实验报告.docx
《偏微分方程数值算法和matlab实验报告.docx》由会员分享,可在线阅读,更多相关《偏微分方程数值算法和matlab实验报告.docx(9页珍藏版)》请在冰点文库上搜索。
偏微分方程数值算法和matlab实验报告
偏微分方程数值实验报告六
实验题目:
用Ritz-Galerkin方法求边值问题
的第n次近似
,基函数为
并用表格列出0.25,0.5,0.75三点处的真解和
时的数值解。
实验算法:
将上述边值问题转化为基于虚功方程的变分问题为:
求
,满足
其中
记
,引入
的n维近似子空间
利用Ritz-Galerkin公式:
,则问题关于
下的近似变分问题解
中的系数
满足
程序代码:
%主函数
pde=modeldata();
I=[0,1];
%积分精度
quadorder=10;
n=[1,2,3,4];
fori=1:
length(n)
uh{i}=Galerkin(pde,I,n(i),quadorder);
end
showsolution(uh{1},'-bo');
holdon
showsolution(uh{2},'-rx');
holdon
showsolution(uh{3},'-.ko');
holdon
showsolution(uh{4},'--k*');
holdoff
title('thesolutionoftheun');
xlabel('x');
ylabel('u');
legend('n=1','n=2','n=3','n=4');
%计算这三点的数值解
x=[1/4,1/2,3/4];
un=zeros(length(n),3);
fori=1:
length(n)
[v,]=basis(x,n(i));
un(i,:
)=(v'*uh{i})';
end
formatshorte
disp('un');
disp(un);
%存储数据
functionpde=modeldata()
pde=struct('f',@f);
functionz=f(x)
z=x.*x;
end
end
%Galerkin方法
functionuh=Galerkin(pde,I,n,quadorder)
h=I
(2)-I
(1);
[lambda,weight]=quadpts1d(I,quadorder);
nQuad=length(weight);
A=zeros(n,n);
b=zeros(n,1);
forq=1:
nQuad
gx=lambda(q);
w=weight(q);
x=[0.25;0.5;0.75];
[phi,gradPhi]=basis(gx,n);
A=A+(-gradPhi*gradPhi'+phi*phi')*w;
b=b+pde.f(gx)*phi*w;
end
A=h*A;
b=h*b;
uh=A\b;
end
%构造基函数
function[phi,gradPhi]=basis(x,n)
m=length(x);
w=sin(pi*x);
v=ones(n,m);
v(2:
end,:
)=bsxfun(@times,v(2:
end,:
),x);
v=cumprod(v,1);
phi=bsxfun(@times,v,w);
gw=pi*cos(pi*x);
gv=[zeros(1,m);v(1:
end-1,:
)];
gv(3:
end,:
)=bsxfun(@times,(2:
n-1)',gv(3:
end,:
));
gradPhi=bsxfun(@times,v,gw)+bsxfun(@times,gv,w);
end
%数值解的图像
functionshowsolution(u,varargin)
x=0:
0.01:
1;
n=length(u);
[v,]=basis(x,n);
y=v'*u;
plot(x,y,varargin{:
});
%varargin:
aninputvariableinafunction
end
%系数矩阵A
function[lambda,weight]=quadpts1d(I,quadorder)
numPts=ceil((quadorder+1)/2);
ifnumPts>10
numPts=10;
end
switchnumPts
case1
A=[02.0000000000000000000000000];
case2
A=[0.57735026918962576450914881.0000000000000000000000000
-0.57735026918962576450914881.0000000000000000000000000];
case3
A=[00.8888888888888888888888889
0.77459666924148337703585310.5555555555555555555555556
-0.77459666924148337703585310.5555555555555555555555556];
case4
A=[0.33998104358485626480266580.6521451548625461426269361
0.86113631159405257522394650.3478548451374538573730639
-0.33998104358485626480266580.6521451548625461426269361
-0.86113631159405257522394650.3478548451374538573730639];
case5
A=[00.5688888888888888888888889
0.538469*********09103631440.4786286704993664680412915
0.90617984593866399279762690.2369268850561890875142640
-0.538469*********09103631440.4786286704993664680412915
-0.90617984593866399279762690.2369268850561890875142640];
case6
A=[0.23861918608319690863050170.4679139345726910473898703
0.66120938646626451366139960.3607615730481386075698335
0.93246951420315202781230160.1713244923791703450402961
-0.23861918608319690863050170.4679139345726910473898703
-0.66120938646626451366139960.3607615730481386075698335
-0.93246951420315202781230160.1713244923791703450402961];
case7
A=[00.4179591836734693877551020
0.40584515137739716690660640.3818300505051189449503698
0.74153118559939443986386480.2797053914892766679014678
0.94910791234275852452618970.1294849661688696932706114
-0.40584515137739716690660640.3818300505051189449503698
-0.74153118559939443986386480.2797053914892766679014678
-0.94910791234275852452618970.1294849661688696932706114];
case8
A=[0.183********564980493947610.3626837833783619829651504
0.52553240991632898581773900.3137066458778872873379622
0.79666647741362673959155390.2223810344533744705443560
0.96028985649753623168356090.1012285362903762591525314
-0.183********564980493947610.3626837833783619829651504
-0.52553240991632898581773900.3137066458778872873379622
-0.79666647741362673959155390.2223810344533744705443560
-0.96028985649753623168356090.1012285362903762591525314];
case9
A=[00.3302393550012597631645251
0.32425342340380892903853800.3123470770400028400686304
0.61337143270059039730870200.2606106964029354623187429
0.83603110732663579429942980.180********48574040584720
0.96816023950762608983557620.0812743883615744119718922
-0.32425342340380892903853800.3123470770400028400686304
-0.61337143270059039730870200.2606106964029354623187429
-0.83603110732663579429942980.180********48574040584720
-0.96816023950762608983557620.0812743883615744119718922];
case10
A=[0.14887433898163121088482600.2955242247147528701738930
0.43339539412924719079926590.2692667193099963550912269
0.67940956829902440623432740.2190863625159820439955349
0.86506336668898451073209670.1494513491505805931457763
0.97390652851717172007796400.0666713443086881375935688
-0.14887433898163121088482600.2955242247147528701738930
-0.43339539412924719079926590.2692667193099963550912269
-0.67940956829902440623432740.2190863625159820439955349
-0.86506336668898451073209670.1494513491505805931457763
-0.97390652851717172007796400.0666713443086881375935688];
end
lambda=(I
(2)+I
(1)+(I
(2)-I
(1))*A(:
1))/2;
weight=A(:
2)*(I
(2)-I
(1))/2;
数值解:
-3.0184e-02-4.2686e-02-3.0184e-02
-2.1919e-02-4.2686e-02-3.8448e-02
-2.3232e-02-4.0652e-02-3.9761e-02
-2.3472e-02-4.0652e-02-3.9521e-02
图像:
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 微分方程 数值 算法 matlab 实验 报告