《2022年微分方程一类边值问题的数值求解》由会员分享，可在线阅读，更多相关《2022年微分方程一类边值问题的数值求解（8页珍藏版）》请在金锄头文库上搜索。
1、学习必备欢迎下载（偏）微分方程一类边值问题的数值求解本文介绍了椭圆型 （偏）微分方程一类边值问题的数值求解程序（笔者自编）及其使用方法。程序基于差分原理，将连续的（偏）微分方程在节点上离散化，最终化为线性代数方程组。 对于原理方面的问题很多微分方程数值解的参考书中都有详尽的描述， 但一般缺少具体的实现程序。 笔者认为， 对这类数值方法是否理解并能运用的关键之一还是在于是否能写出用于解决问题的程序。理论基础固然必不可少，程序确往往是我们解决问题的敲门砖。一维椭圆型方程可表示如下：,0)(,0)()(baxxqxpfqudxdurdxdupdxdLu其中 L 表示微分算子，很明显这是一个线性算子。
2、这里要求p，q 大于零是为了保证最终得到的线性方程组有唯一的非零解，但事实上不满足这个条件可能也是有解的， 这涉及到微分方程解的存在性也确定性问题，读者若有兴趣可参考相关书籍。更具体的，我们可以举出一个一维椭圆型方程的例子来：例（1） ：)cos(10)exp(222xuxdxduxdxud即：)cos(10),exp()(,)(,1)(2xfxxqxxrxp3)3(1)1 (uu利用文中提及的程序，可以将以上问题表述为：syms x %定义符号变量 xp=(x) 1;r=(x) x.2;q=(x) exp(x);f=(x) 10*cos(x);Odbvp(p,r,q,f,1,3,1,3);
3、精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 1 页，共 8 页学习必备欢迎下载11.21.41.61.822.22.42.62.83-0.500.511.522.532Odbvp u=f(x)xu哈哈，一条非常漂亮的曲线。若不满足 p,q0，我们也举一例：例（2）syms xp=(x) 10*cos(x);r=(x) x.2;q=(x) 10*sin(x);f=(x) 10*cos(x);Odbvp(p,r,q,f,1,3,1,3) 11.21.41.61.822.22.42.62.83-500501001502002503003502Odbv
4、p u=f(x)xu精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 2 页，共 8 页学习必备欢迎下载可以发现，程序仍能求解， 但结果的光滑性不好。 这可能与方程本身的刚性有关，也可能和程序所适用的范围有关。事实上，在求解线性方程组时用到了matlab中的反除原理（ Ab） ，这个方法本质上是基于最小二乘原理的拟合方法，因此也可以用来求解线性方程组（即便方程是矛盾的，奇异的）。因此，若用上笔者所给的程序计算， 总是能算出一个结果， 但结果的正确性与合理性要依靠使用者对问题的判断，在满足p,q0的前提下，解是没有问题的。更进一步的，对于二维椭圆型方程
5、可表示如下：,0)(,0)()()(baxxqxpfquyupyxupxLu二维边值的一类边界条件又可根据边界的形状化分为多个类型， 为简单起见此处介绍矩形边界一类边值问题的求解。如边界固定的混凝土板受均布荷载作用时的微分方程为：例（3）fyuxu2222clc ;clear;syms xyp=(x,y) 1;q=(x,y) 0;f=(x,y) -0.3;BF=(x,y) 0;Pdbvp(p,q,f,-3,3,-3,3,BF,40,40) ;计算所得板各点挠度如下图：精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 3 页，共 8 页学习必备欢迎下载
6、-3-2-10123-202-0.8-0.6-0.4-0.20x2Pdbvpyu-0.7-0.6-0.5-0.4-0.3-0.2-0.10再看一个复杂点的问题，二维问题的求解相对较难，但结果往往非常有趣。我们可以得到一些非常有意思的图像。在看到一个方程时， 它只是一段数学表达式，它一般有明确的物理意义， 但对于我们来说它始终如黑匣子一般，大多时候我们是不知其庐山真面目的。 笔者常常有一种渴望， 渴望看看它们的真实样子， ，这也正是促使笔者做这个工作的原因之一，如：例（4）clc ;clear;syms xyp=(x,y) abs(sin(x)+exp(y)*0.002);q=(x,y) y.2
7、-x.2;f=(x,y) -sin(x)+cos(x);BF=(x,y) 0;Pdbvp(p,q,f,-3,3,-3,3,BF,40,40) ;精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 4 页，共 8 页学习必备欢迎下载-202-202-6-4-2024x2Pdbvpyu-5-4-3-2-1012345又一张非常漂亮的曲面！上述计算程序的专用性非常强， 使问题的解决更加简洁， 往往只需要几行程序就能解决这类边值问题，使用起来直观易懂，非常适于初学者使用。当然，专用性强就必然导致通用性差， 上述椭圆型方程只是众多微分方程中的一个非常小的部分，而
8、所使用的差分方法也是精度较低的一种，对上述范围以外的问题是无能为力的。由于笔者偷懒，对程序中产生的稀疏矩阵未加处理，这对于求解是不利的。在节点划分增加时， 如若不处理稀疏矩阵求解可能根本进行不下去。对于大规模的线性方程组，反除法也是不适用的，迭代法、共轭梯度法往往更加高效。笔者写这个也是突然兴起， 对这类问题进行了一个小的尝试，问题不少当然趣味也不少，只希望对大家有一定的帮助 附程序一维：function u,xx=Odbvp(p,r,q,f,a,b,alpha,beta,n) % Directory%Ordinary differential problem based on the fir
9、st boundray condition(can be used one or second order odbvp problem).%p(x),r(x),q(x),f(x) are functions based on the practical problems. %a,bis the interval which need to be caculated.%alpha,beta is the boundray value.%n is the number of the nodes in the inteval a,b% Grid if nargin9精选学习资料 - - - - -
10、- - - - 名师归纳总结 - - - - - - -第 5 页，共 8 页学习必备欢迎下载 n=50;endD=zeros(n+1,n+1);G=zeros(n+1,1);u=G;h=(b-a)/n;syms xX=(x) a+(b-a)/(n+1)*x;% Caculate The Matix G&Dfor i=1:n+1if i=1 G(1)=alpha; D(i,i)=1;endif i=n+1 G(i)=beta; D(i,i)=1;endif i=1&i=n+1 A(i)=2*feval(p,X(i-1/2)./h+feval(r,X(i); C(i)=2*feval(p,X(i
11、+1/2)./h-feval(r,X(i); B(i)=2*(feval(p,X(i-1/2)+feval(p,X(i+1/2)./h+2*h*feval(q,X(i); G(i)=2*h*feval(f,X(i); D(i,i-1)=-A(i); D(i,i)=B(i); D(i,i+1)=-C(i);endend% Solve The Linear Equationsu=DG;% Plotxx=a:h:b;plot(xx,u);grid on;title(2Odbvp u=f(x);xlabel(x);ylabel(u);end二维：function U,Xx,Yy=Pdbvp(p,q,f
12、,ax,bx,ay,by,BF,n1,n2) % Directory%Partial differential problem with the first bundry condition.%p(x,y),q(x,y),f(x,y) are functions based on practical problems;%ax,bx*ay,by is area need to be caclulated.精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 6 页，共 8 页学习必备欢迎下载%BF(x,y)is boundry condition func
13、tion;%n1,n2 is the mesh desity in two dimensions.% Grid if nargin10 n1=30; n2=30;endsyms xyX=(x) ax+(bx-ax)/(n1+1)*x;Y=(y) ay+(by-ay)/(n2+1)*y;h1=(bx-ax)/n2;h2=(by-ay)/n1;H=zeros(n1+1)*(n2+1),(n1+1)*(n2+1);F=zeros(n1+1)*(n2+1),1);% Caculate The Matix H&FNd=0;for j=1:n2+1for i=1:n1+1 Nd=Nd+1;if i=1|i
14、=n1+1|j=1|j=n2+1 H(Nd,Nd)=1; F(Nd)=BF(X(i),Y(j);else a(1)=feval(p,X(i+1/2),Y(j)./h1.2; a(2)=feval(p,X(i),Y(j+1/2)./h2.2; a(3)=feval(p,X(i-1/2),Y(j)./h1.2; a(4)=feval(p,X(i),Y(j-1/2)./h2.2; a0=sum(a)+feval(q,X(i),Y(j); H(Nd,Nd-n1-1)=-a(4); H(Nd,Nd-1)=-a(3); H(Nd,Nd)=a0; H(Nd,Nd+1)=-a(1); H(Nd,Nd+n1+1
15、)=-a(2); F(Nd)=feval(f,X(i),Y(i);endendend% Solve The Linear Equationsu=HF;% Data SortingXx,Yy=meshgrid(ax:h1:bx,ay:h2:by);Nd=0;for j=1:n2+1for i=1:n1+1精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 7 页，共 8 页学习必备欢迎下载 Nd=Nd+1; U(i,j)=u(Nd);endendH;% Surfsurfc(Xx,Yy,U);shading interp%meshc(Xx,Yy,U);colormap jetxlabel(x);ylabel(y);zlabel(u);axis equalcolorbar;title(2Pdbvp)%hold on%C0,h = contour3(Xx,Yy,U); %clabel(C0,h);%contour(Xx,Xy,U,0);end精选学习资料 - - - - - - - - - 名师归纳总结 - - - - - - -第 8 页，共 8 页
