当前位置：文档之家 › 二阶椭圆偏微分方程实例求解(附matlab代码)

二阶椭圆偏微分方程实例求解(附matlab代码)

  • 格式：doc
  • 大小：254.40 KB
  • 文档页数：13

《微分方程数值解法》期中作业实验报告二阶椭圆偏微分方程第一边值问题姓名：学号：班级：2013年11月19日二阶椭圆偏微分方程第一边值问题摘要对于解二阶椭圆偏微分方程第一边值问题，课本上已经给出了相应的差分方程。

而留给我的难题就是把差分方程组表示成系数矩阵的形式，以及对系数进行赋值。

解决完这个问题之后，我在利用matlab 解线性方程组时，又出现“out of memory ”的问题。

因为99*99阶的矩阵太大，超出了分配给matlab 的使用内存。

退而求其次，当n=10，h=1/10或n=70，h=1/70时，我都得出了很好的计算结果。

然而在解线性方程组时，无论是LU 分解法或高斯消去法，还是gauseidel 迭代法，都能达到很高的精度。

关键字：二阶椭圆偏微分方程 差分方程 out of memory LU 分解 高斯消去法 gauseidel 迭代法一、题目重述解微分方程：()()2222((,))((,))()(,)()(,)(,)1y x x x y y x y yxxyxye u x y e u x y x y u x y x y u x y u x y y e x e e y x e--+++-+=-++++已知边界:(0,)1,(1,),(,0)1,(,1)y x u y u y e u x u x e ====求数值解, 把区域[0,1][0,1]G =?分成121/100,1/100h h ==，n =100 注：老师你给的题F 好像写错了，应该把22x y y e x e +改成22y x y e x e +。

二、问题分析与模型建立2.1微分方程上的符号说明()()22221y x xy xy y e x e e y x e -++++2.2课本上差分方程的缺陷课本上的差分方程为：举一个例子：当i=2,j=3时，；当i=3,j=3时，。

但，是，显然这两个不是同一个数，其大小也不相等。

2.3差分方程的重新定义因此，为了避免2.2中赋值重复而产生的错误，我在利用matlab编程时，对这些系数变量进行了重新定义：，，，，2.4模型建立这里的模型建立就是把差分方程组改写成系数矩阵的形式。

经过研究，我觉得写成如下的系数矩阵不仅看起来简单明了，而且在matlab编程时比较方便。

系数矩阵为：Pu=f其中P是（）阶方阵，具体如下：而u是（）维的列向量，具体如下：而f是（）维的列向量，具体如下：三、求解过程3.1对系数矩阵的分析对上述模型的求解就是对线性方程组的求解。

通过观察，我发现P是一个对角占优的矩阵，这不仅确定了解的唯一性，还保证了迭代法的收敛性。

此外，还可以确定进行LU分解，若使用高斯消去法还可以省去选主元的工作。

3.2matlab编程因为矩阵阶数过大，所以此题的编程难点为构造系数矩阵，即对线性方程组的赋值。

我采用的方法是分块赋值。

对于P的赋值，过程如下：第一步：第二步：第三步：P=BCD-diag(G,99)-diag(K,99).P和 f的具体构造见附录6.1主代码3.3编程求解过程中的问题3.3.1问题产生当按照老师要求，n=100，h=1/100时，运行编好的matlab程序时，会出现如图3.1的错误提示。

图3.13.3.2问题分析在matlab的命令窗口输入“memory”，出现如图3.2的内存使用情况，可以得出：Memory used by MATLAB: 454 MB (4.760e+008 bytes)。

，若不用稀疏矩阵定义P，经过粗略计算，我发现矩阵P就要占800MB左右的内存，加上其他数据，内存消耗至少在1G以上。

但是我电脑上分配给matlab的内存只有：454 MB，即使在关闭杀毒软件等大部分应用程序后，分配给matlab的内存也刚够1G。

图3.23.3.3问题解决经过上网查找资料后，我找到了如下几个解决方法。

1）尽量早的分配大的矩阵变量2）尽量避免产生大的瞬时变量，当它们不用的时候应该及时clear。

3）尽量的重复使用变量（跟不用的clear掉一个意思）4）将矩阵转化成稀疏形式5）使用pack命令6）如果可行的话，将一个大的矩阵划分为几个小的矩阵，这样每一次使用的内存减少。

7）增大内存针对本题，我觉得比较理想的解决方法是采用稀疏矩阵的方式定义P。

这样可以有效的减小P的内存消耗。

但是考虑到老师的这次期中作业主要是考察我们对二阶椭圆偏微分方程的理解与实例操作，而不是旨在考察我们的matlab编程能力。

因此我在此，略作偷懒，把n改成10或70（75以上内存就不够用了），适当的降低精度来得到结果。

四、计算结果4.1当n=10,h=1/10时的结果取n=10,h=1/10时，matlab运行的部分结果如图4.1。

表4.2为LU分解法和高斯消去法的部分结果（这两个直接法结果完全一样），表4.3为迭代法的部分结果。

图4.14.2当n=70,h=1/70时的结果取n=70,h=1/70时，matlab运行的部分结果如图4.4（LU分解法）。

计算时间为17多分钟（1043秒），误差至少在小数点后9位。

图4.4五、结论分析由于本人的电脑比较破旧，内存不是很大，再加上没有采取稀疏矩阵的存储方式，难以达到n=100，h=1/100的计算要求。

但改为n=10，h=1/10或n=70，h=1/70后，计算精度也十分理想。

尤其是后者，其误差至少在小数点后9位，在比较使用哪种方法解线性方程组时，直接法中的LU分解法和高斯消去法的计算结果是完全相等的。

而gauseidel迭代法计算结果按个人设定的计算精度而定。

对于这种大型线性方程组来说，迭代法从计算速度和计算机存储方面来看具有超过直接法的决定性优点。

对于稀疏方程组来说，迭代法十分有效。

从本题的实际情况来看，当n=70时，LU分解的直接法的计算时间为17分钟左右，而gauseidel迭代法的计算时间为20秒左右，这与以上分析的情况一致。

但是当n 越来越大时（从n=10起），固定迭代步数的迭代法解的精度越来越差，为了达到与直接法相同的精度，必须提高迭代步数。

然而这又会加大计算量，使计算速度变慢（见表5.1）仅从本题计算结果来看（n=10时）：当x，y靠近0时，直接法的数值解更准确，而当x，y靠近1时，迭代法的数值解更准确。

这与gauseidel迭代法的算法特点相符合。

因为gauseidel迭代法后面的解在迭代时要用到前面的解，所以x，y靠近1的后面的解更为精确。

六、附录6.1主代码主代码中解决了对系数矩阵的赋值，即写出了线性方程组。

在解线性方程组时可以选用LU分解法、高斯消去法和迭代法中的任意一种。

function n=Middle_Term_Bymyself(t)%t为区间划分数tic;%定义函数及初始化基本变量FA=@(x,y)exp(y);FB=@(x,y)exp(x);FC=@(x,y)x+y;FD=@(x,y)x-y;FE=@(x,y)1;FF=@(x,y)(y^2+x^2+1)*exp(x*y)-(y^2*exp(y)+x^2*exp(x))*exp(x*y);FU=@(x,y)exp(x*y);%真实值n=t;%区间划分为n*nh=1/n;%h为单位区间长度A=zeros(2*n-1);B=zeros(2*n-1);C=zeros(n-1);D=zeros(n-1);E=zeros(n-1);F=zeros(n-1);U=zeros(n-1);%真实值的矩阵表示u=zeros((n-1)*(n-1),1);%真实值的数列表示error=zeros((n-1)*(n-1),1);%误差P=zeros((n-1)*(n-1));%最终要解的方程组的系数矩阵,即平时的A f=zeros((n-1)*(n-1),1); %即平时的bff=zeros(n-1); %ff(i,j)=f((i-1）*9+j)BCD=[];%记录三对角部分G=[];%记录上三角的那一斜条K=[];%记录下三角的那一斜条b=zeros(n-1);%表示a(i,j)c=zeros(n-1);%表示a(i,j+1)d=zeros(n-1);%表示a(i,j-1)g=zeros(n-1);%表示a(i+1,j)k=zeros(n-1);%表示a(i-1,j)%对函数进行赋值x=zeros(2*n-1,1);y=zeros(2*n-1,1);for i=1:2*n-1 %使A.B下标i-1/2变为2i-1for j=1:2*n-1x(i)=i*h/2;y(j)=j*h/2;A(i,j)=FA(x(i),y(j));B(i,j)=FB(x(i),y(j));endendx=zeros(n-1,1);y=zeros(n-1,1);for i=1:n-1for j=1:n-1x(i)=i*h;y(j)=j*h;C(i,j)=FC(x(i),y(j));D(i,j)=FD(x(i),y(j));E(i,j)=1;F(i,j)=FF(x(i),y(j));U(i,j)=FU(x(i),y(j));endend%对最终要解的方程组的系数矩阵进行赋值for i=1:n-1bcd=zeros(n-1);bb=[];cc=[];dd=[];gg=[];kk=[];for j=1:n-1b(i,j)=(A(2*i+1,2*j)+A(2*i-1,2*j)+B(2*i,2*j+1)+B(2*i,2*j-1))/h^2+E(i,j);c(i,j)=(B(2*i,2*j+1)-h*D(i,j)/2)/h^2;d(i,j)=(B(2*i,2*j-1)+h*D(i,j)/2)/h^2;g(i,j)=(A(2*i+1,2*j)-h*C(i,j)/2)/h^2;k(i,j)=(A(2*i-1,2*j)+h*C(i,j)/2)/h^2;bb=[bb b(i,j)];if j<=n-2cc=[cc c(i,j)];endif j>=2dd=[dd d(i,j)];endgg=[gg g(i,j)];kk=[kk k(i,j)];%给f赋值if i==1ff(i,j)=F(i,j)+k(i,j)*1;%边值为1elseif i==n-1ff(i,j)=F(i,j)+g(i,j)*A(i,2*j);%A中i取值无所谓，不影响elseff(i,j)=F(i,j);endif j==1ff(i,j)=ff(i,j)+d(i,j)*1;%边值为1elseif j==n-1ff(i,j)=ff(i,j)+c(i,j)*B(2*i,j);endf((i-1)*(n-1)+j,1) = ff(i,j);%你应该懂的，坐标变换endbcd=diag(bb)-diag(cc,1)-diag(dd,-1);BCD=blkdiag(BCD,bcd);if i<=n-2G=[G gg];endif i>=2K=[K kk];endendP=BCD-diag(G,n-1)-diag(K,-n+1);%BCD=BCD-diag(G,n-1)-diag(K,-n+1);x=Doolittle(BCD,f);这样是不是可以减少点内存消耗x=Doolittle(P,f);%LU分解解方程组，这里也可以输入其他解方程组的方法%x=GaussXQByOrder(P,f) %高斯消去法%x0=ones((n-1)^2,1);x=gauseidel(P,f,x0) %gauseidel迭代法for i=1:n-1for j=1:n-1u((i-1)*(n-1)+j)=U(i,j);endenderror=abs(x-u);result=[x u error];time = toc;disp('计算时间为：');disp(time);disp('---------------------------------------------------------------');format long;disp('计算结果为：');disp(' 数值解真实值误差');disp(result);6.2LU分解解线性方程组function [x,L,U]= Doolittle (A,b)N = size(A);n = N(1);L = eye(n,n); %L的对角元素为1U = zeros(n,n);U(1,1:n) = A(1,1:n); %U的第一行L(1:n,1) = A(1:n,1)/U(1,1); %L的第一列for k=2:nfor i=k:nU(k,i) = A(k,i)-L(k,1:(k-1))*U(1:(k-1),i);%U的第k行endfor j=(k+1):nL(j,k) = (A(j,k)-L(j,1:(k-1))*U(1:(k-1),k))/U(k,k);%L的第k列endendy = SolveDownTriangle(L,b);x = SolveUpTriangle(U,y); %求解方程6.3高斯消去法function [x,XA]=GaussXQByOrder(A,b)N = size(A);n = N(1);for i=1:(n-1)for j=(i+1):nif(A(i,i)==0)disp('对角元素为0！'); %防止对角元素为0return;endl = A(j,i);m = A(i,i);A(j,1:n)=A(j,1:n)-l*A(i,1:n)/m; %消元方程b(j)=b(j)-l*b(i)/m;endendx=SolveUpTriangle(A,b); %通用的求上三角系数矩阵线性方程组的函数XA = A; %消元后的系数矩阵6.4上三角解线性方程组（LU分解法、高斯消去法要用到这个算法）function x=SolveUpTriangle(A,b)N=size(A);n=N(1);for i=n:-1:1if(i<n)s=A(i,(i+1):n)*x((i+1):n,1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end6.5下三角解线性方程组（LU分解法要用到这个算法）function x=SolveDownTriangle(A,b)N=size(A);n=N(1);for i=1:nif(i>1)s=A(i,1:(i-1))*x(1:(i-1),1);elses=0;endx(i,1)=(b(i)-s)/A(i,i);end6.6gauseidel迭代法if nargin==3eps= 1.0e-6;%可以安自己的精度要求改变M = 10000;elseif nargin == 4M = 10000;elseif nargin<3errorreturn;endD=diag(diag(A)); %求A的对角矩阵L=-tril(A,-1); %求A的下三角阵U=-triu(A,1); %求A的上三角阵G=(D-L)\U;f=(D-L)\b;x=G*x0+f;n=1; %迭代次数while norm(x-x0)>=epsx0=x;x=G*x0+f;n=n+1;if(n>=M)disp('Warning: 迭代次数太多，可能不收敛！'); return;endenddisp('迭代次数为：');disp(n);。

合集下载

偏微分方程Matlab数值解法(补充4)

偏微分方程Matlab数值解法(补充4)

偏微分方程Matlab 数值解法（补充4）Matlab 可以求解一般的偏微分方程，也可以利用偏微分方程工具箱中给出的函数求解一些偏微分方程。

1 偏微分方程组求解Matlab 语言提供了pdepe()函数，可以直接求解偏微分方程(,,,)[(,,,)](,,,)m mu u u u C x t u x x f x t u s x t u x t x x x-∂∂∂∂∂=+∂∂∂∂∂ （4.1）这样，偏微分方程可以编写为以下函数的描述，其入口为[,,](,,,)x c f s pdefun x t u u =其中：pdefun 为函数名。

由给定输入变量可计算出,,c f s 这三个函数。

边界条件可以用下面的函数描述(,,)(,,).*(,,,)0up x t u q x t u f x t u x∂+=∂ （4.2） 这样的边值函数可以写为Matlab 函数[,,,](,,,)a a b b x p q p q pdebc x t u u =初始条件数学描述为00(,)u x t u =，编写一个简单的函数即可0()u pdeic x =还可以选择x 和t 的向量，再加上描述这些函数，就可以用pdepe （）函数求解次偏微分方程，需要用如下格式求解(,@,@,@,,)sol pdepe m pdefun pdeic pdebc x t =【例1】 试求下列偏微分方程2111222221220.024()0.17()u u F u u t x u u F u u tx ⎧∂∂=--⎪⎪∂∂⎨∂∂⎪=+-⎪∂∂⎩ 其中： 5.7311.46()xx F x e e -=-，且满足初始条件1(,0)1u x =，2(,1)0u x =及边界条件1221(0,)0,(0,)0,(1,)1,(1,)0u u t u t u t t x x∂∂====∂∂解：对照给出的偏微分方程和（4.1），可将原方程改写为111222120.024/()1.*10.17/()u u x F u u u u x F u u t x ∂∂--⎡⎤⎡⎤⎡⎤⎡⎤∂∂=+⎢⎥⎢⎥⎢⎥⎢⎥∂∂-∂∂⎣⎦⎣⎦⎣⎦⎣⎦可知0m =，且1122120.024/()1,,10.17/()u x F u u c f s u x F u u ∂∂--⎡⎤⎡⎤⎡⎤===⎢⎥⎢⎥⎢⎥∂∂-⎣⎦⎣⎦⎣⎦编写下面的Matlab 函数function [c,f,s]=c7mpde(x,t,u,du)c=[1;1];y=u(1)-u(2);F=exp(5.73*y)-exp(-11.46*y);s=F*[-1;1]; f=[0.024*du(1);0.17*du(2)];套用（4.2）中的边界条件，可以写出如下的边值方程左边界2010.*00f u ⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦，右边界1100.*100u f -⎡⎤⎡⎤⎡⎤+=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦⎣⎦ 编写下面的Matlab 函数function [pa,qa,pb,qb]=c7mpbc(xa,ua,xb,ub,t) pa=[0;ua(2)];qa=[1;0];pb=[ub(1)-1;0];qb=[0,1]; 另外，描述初值的函数function u0=c7mpic(x) u0=[1;0];有了这三个函数，选定x 和t 向量，则可以由下面的程序直接求此微分方程，得出解1u 和2u 。

matlab 求解椭圆方程

matlab 求解椭圆方程

matlab 求解椭圆方程椭圆方程是数学中的一个重要概念，它描述了平面上的一种几何形状。

在数学中，椭圆是一个集合，它包含了与两个给定点P和Q的距离之和等于常数2a的所有点的集合。

这个常数2a被称为椭圆的长轴，而P和Q被称为椭圆的焦点。

在实际应用中，解椭圆方程是非常有用的。

例如，在物理学中，椭圆方程可以用来描述行星绕太阳的轨道。

在工程学中，椭圆方程可以用来设计椭圆形的天线。

在经济学中，椭圆方程可以用来描述收入分配的不平等程度。

为了解椭圆方程，我们可以使用MATLAB这个功能强大的数学软件。

MATLAB提供了一个名为"ellipticE"的函数，可以用来计算椭圆的椭圆积分。

椭圆积分是计算椭圆弧长、面积和体积等参数的重要工具。

我们需要定义椭圆的参数。

这包括椭圆的长轴a和短轴b。

然后，我们可以使用"ellipticE"函数来计算椭圆的椭圆积分。

在MATLAB 中，我们可以这样写代码：```a = 3; % 椭圆的长轴b = 2; % 椭圆的短轴% 计算椭圆的椭圆积分E = ellipticE(pi/4, a/b);```在这个例子中，我们计算了一个长轴为3，短轴为2的椭圆的椭圆积分。

我们使用了"pi/4"作为参数，这代表了椭圆的弧度。

计算得到的结果存储在变量E中。

除了椭圆积分，MATLAB还提供了其他与椭圆相关的函数，如椭圆的周长、面积和体积等。

我们可以根据具体的需求选择适当的函数来解决问题。

椭圆方程是数学中的一个重要概念，它在实际应用中有着广泛的应用。

通过使用MATLAB提供的函数，我们可以方便地解决椭圆方程，并计算椭圆的各种参数。

这为我们在科学研究和工程设计中提供了有力的工具和方法。

通过深入学习和理解椭圆方程的性质和特点，我们可以更好地应用它们解决实际问题。

偏微分方程中的椭圆方程的matlab解法

偏微分方程中的椭圆方程的matlab解法

偏微分方程中的椭圆方程的matlab解法
摘要：对偏微分方程中的椭圆方程的matlab解法进行研究，给出几个椭圆型方程的例子，然后用matlab软件进行求解。

分析结果，比较结果。

关键词：偏微分方程；椭圆方程；matlab
一问题的提出
椭圆型方程已经被许多学者所关注，关于其的边值得问题的研究有很多的工作，已经有了不少的结果，我们看一下。

在xy平面的某区域Ω中，未知函数u满足Laplace方程，将边界分成若干弧段，要求u在每一弧段上满足不同类型的边界条件。

讨论此类定解问题的差分模拟。

例如，求解如下定解问题：
类似于此类的问题有很多很多种，我们在这里考虑一些比较简单的问题。

二记号及引理
引理 1 若令 c ( x ) ∈ C ( Ω), c ( x ) ≥ 0 且 c ( x ) 不恒为零，β ( x ) ≡ 0 ，x ∈∂Ω, 假使 u ∈ C (Ω) I C ( Ω ) 且满足下述不等式
Li u + c( x)u ≥ 0 x ∈Ω,
B u ( x) ≥ 0 x ∈∂Ω,
则u ( x) ≥ 0 ，∀x ∈Ω .
进而，若Li u + c( x)u不恒等于0，
或者x ∈Ω，Bi u ( x)不恒等于0， x ∈∂Ω ,
则u ( x, t ) > 0 ( x, t ) ∈Ω×[0, ω ].
引理 2 α(i = 1,2,L, n)，如果 f i 满足 ( F1 ) ，那么 H i 和 hi 也满足 ( F1 ) .
三主要结果及证明。

matlab解二阶偏微分方程

matlab解二阶偏微分方程

<div align="center">一、Matlab求解二阶偏微分方程(ODE)的基本步骤</div>1. 数学模型：首先要确定求解的方程是哪一类的偏微分方程(ODE)，然后建立其对应的数学模型，使其符合这类微分方程的形式;2. 确定边界条件：确定迭代范围$[a,b]$,边界条件函数 $y(a)=\alpha$ 、$y(b)=\beta$;3. 写出Matlab程序：在该类ODE中，通常会有某一种常用的数值求解方法，一般使用微分方程求解器(ODE)，如ode45等;4. 获得实际结果：开始编写Matlab程序，完成参差和参数的输入以后，可以运行Matlab程序，然后求得结果，再用图像表示出来。

<div align="center">二、具体求解</div>$$\frac{d^2y}{dx^2}+y=6sin(2x)$$微分方程为二阶常微分方程，求解条件如下：$[a,b]=[0,\pi], y(0)=1,y(\pi)=3.$（1）Matlab函数表达式首先建立与二阶非齐次线性常微分方程相符合的数学模型，其Matlab函数表达式为$$ f(x,y,y')=\frac{dy}{dx}-y'-6sin2x $$其中，$y=y(x)$；（2）函数程序在Matlab中，定义函数程序 $myode.m$ ，此程序返回右端函数 $f(x,y,y')$ 的值表达式，程序内容如下。

```MATLAB% 右端函数程序function dy=myode(x,y)dy=[y(2);-y(2)-6*sin(2*x)];end```（3）调用Matlab函数olvede45调用Matlab函数 solvede45 求解二阶ODN，程序内容如下：```MATLAB% 主程序求解% maxstep表示分裂的步长大小% Tolerence表示误差，控制求解精度Maxstep=0.25;Tolerence=1e-4;a=0;b=pi;y0=[1;0];[x,y] = ode45('myode',[a,b],y0,options);```（4）结果展示输入参数之后，运行Matlab程序，得到如下图：![](../images/matlab_2_diff.png)此图为$y(x)$随$x$变化的曲线，可以看出，二阶偏微分的求解结果满足了边界条件，即$y(0)=1,y(\pi)=3$ ，如图中红色圆点所示。