偏微分方程数值算法和matlab实验报告

偏微分方程数值实验报告六

实验题目：

用Ritz-Galerkin 方法求边值问题

2

(0)0,(1)1

u u x u u ''-+===⎧⎨⎩01x << 的第n 次近似()n u x ，基函数为(),i 1,2,.()..x n i sin i x πφ==并用表格列出0.25，0.5,0.75三点处的真解和1,2,3,4n =时的数值解。 实验算法：

将上述边值问题转化为基于虚功方程的变分问题为： 求

10()u H I ∈，满足1

0(,)(,(v H ,))a u v x x I ∀∈=-

其中

11

00

(u''v uv)(,)(,)dx (u'v'uv)dx a u v Lu v +=-+==⎰⎰ 记(x)sin()w x π=，引入10()H I 的n 维近似子空间12{,},,n n U φφφ=⋯(),i 1,2,,i sin i x n πφ==⋯利用Ritz-Galerkin 公式：1()(f,),j 1,,,2,n j i j

j a n φφφ==⋯=∑，则问题关于n U 下的近似变分问题解

1()()n i i n i c u x x φ==∑中的系数12(,,,)T n

n c c c c =⋯∈满足12()(,,)n j i j j a x φφφ==∑

程序代码：

%主函数

pde=modeldata();

I=[0,1];

%积分精度

quadorder=10;

n=[1,2,3,4];

for i=1:length(n)

uh{i}=Galerkin(pde,I,n(i),quadorder); end

showsolution(uh{1},'-bo');

hold on

showsolution(uh{2},'-rx');

hold on

showsolution(uh{3},'-.ko');

hold on

showsolution(uh{4},'--k*');

hold off

title('the solution of the un');

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);

for i=1:length(n)

[v, ]=basis(x,n(i));

un(i,:)=(v'*uh{i})';

end

format shorte

disp('un');

disp(un);

%存储数据

function pde=modeldata()

pde=struct('f',@f);

function z=f(x)

z=x.*x;

end

end

%Galerkin方法

function uh=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);

for q=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

%数值解的图像

function showsolution(u,varargin)

x=0:0.01:1;

n=length(u);

[v, ]=basis(x,n);

y=v'*u;

plot(x,y,varargin{:});

% varargin: an input variable in a function

end

%系数矩阵A

function [lambda,weight] = quadpts1d(I,quadorder)

numPts = ceil((quadorder+1)/2);

if numPts> 10

numPts = 10;

end

switch numPts

case 1

A = [0 2.0000000000000000000000000];

case 2

A = [0.5773502691896257645091488

1.0000000000000000000000000

-0.5773502691896257645091488

1.0000000000000000000000000];

case 3

A = [ 0

0.8888888888888888888888889

0.7745966692414833770358531

0.5555555555555555555555556