数组与矩阵---未排序数组中累加和小于或等于给定值的最长子数组问题

在给定无序数组arr和整数k的情况下,找到所有子数组中累加和小于或等于k的最长子数组长度。方法一是通过生成累加和数组并使用二分查找,时间复杂度O(NlogN),空间复杂度O(N);方法二是以O(N)的时间复杂度和空间复杂度解决该问题。

摘要生成于 C知道 ,由 DeepSeek-R1 满血版支持, 前往体验 >

【题目】

  给定一个无序数组arr,其中元素可正、可负、可0,给定一个整数k。求arr所有的子数组中累加和小于或等于k的最长子数组长度。
  
  例如:arr = [3, -2, -4, 0, 6],k = -2,相加和小于或等于-2的最长子数组为[3, -2, -4, 0],所以结果返回4。

【基本思路】

  方法一:时间复杂度O(NlogN),空间复杂度O(N)
  
  依次求以数组每个位置结尾的累加和小于或等于 k 的最长子数组长度。如何找呢?首先生成每个位置的累加和数组 sumArr,对于每一个位置 i,要求以该位置结尾的累加和小于或等于k的最长子数组,实际上只要找到 sumArr[0…i] 中第一个大于或等于sum-k的位置,假设位置为 index,那么 arr[index+1…i] 就是以 i 位置结尾的累加和小于或等于 k 的最长子数组。

  那么如何快速的求出 sumArr[0…i] 中第一个大于或等于 sum-k 的位置呢?假设累加和数组sumArr = [0,1, 3, 2, 7, 5]。注意这里数组的第一项0表示没有任何一个数时的累加和,它的存在是为了防止以 0 位置开头的子数组的丢失(这个在未排序数组中累加和为给定值的最长子数组中已经分析过)。对于数组sumArr我们可以将它转变为 [0, 1, 3, 3, 7, 7],为什么呢?因为我们只关心第一个大于 sum-k 的位置,如果累加和为 2 已经大于 sum-k,那么累加和 3 必定也大于 sum-k,所以只要保留一个更大的、出现更早的累加和就可以。转变后的数组很显然是一个有序数组,那么我们就可以使用二分查找在 logN 的时间复杂度内找到第一个大于或等于 sum-k 的位置。如果不转变成有序数组,很显然,我们需要用 N 的复杂度去遍历寻找满足条件的位置。

  具体实现参见如下代码:

#python3.5
def maxLen(arr, k):
    def getLessIndex(h, num, index):
        left = 0
        right = index
        res = 1
        while left <= right:
            mid = (left + right) // 2
            if h[mid] >= num:
                res = mid
                right = mid - 1
            else:
                left = mid + 1
        return res

    if arr == None or len(arr) == 0:
        return 0
    sum = 0
    length = 0
    h = [0 for i in range(len(arr)+1)]
    for i in range(1, len(h)):
        sum += arr[i-1]
        h[i] = max(h[i-1], sum)
    sum = 0
    for i in range(len(arr)):
        sum += arr[i]
        pre = getLessIndex(h, sum-k, i)  #这个位置是h中的下标,注意转换成arr中的下标
        tmpLen = 0 if pre == -1 else i-pre+1
        length = max(length, tmpLen)
    return length

  方法二:时间复杂度O(N),空间复杂度O(N)。

#python3.5
def maxLen2(arr, k):
    if arr == None or len(arr) == 0:
        return 0
    help = [0 for i in range(len(arr))]   #以i位置开头能达到的最小累加和
    help[-1] = arr[-1]
    map = {}   #以位置i开头的最小累加和数组的最右位置
    map[len(arr)-1] = len(arr) - 1
    for i in range(len(arr)-2, -1, -1):
        if help[i+1] <= 0:
            help[i] = help[i+1] + arr[i]
            map[i] = map[i+1]
        else:
            help[i] = arr[i]
            map[i] = i
    res = 0   #全局变量记录最大子数组长度
    end = 0   #子数组右边界的下一个位置
    sum = 0   #子数组累加和
    for i in range(len(arr)):
        #计算以i位置开头满足条件的最大子数组
        while end < len(arr) and sum + help[end] <= k:
            sum += help[end]
            end = map[end] + 1
        #更新最大子数组长度
        res = max(res, end - i)
        #子数组左边界向右移动一位
        sum -= arr[i] if end > i else 0
        #更新end
        end = end if end > i else i + 1
    return res
<think>为了求解给定的二阶常微分方程边值问题,我们将使用有限元方法(FEM)。问题描述如下:-u''+u=(1+π²)sin(πx),在区间(0,1)上边界条件:u(0)=0,u(1)=0已知精确解为u(x)=sin(πx)。有限元方法步骤:1.将区间[0,1]离散为N+1个节点(即N个单元)。2.选择线性基函数(分片线性函数)。3.建立刚度矩阵和载荷向量。4.应用边界条件。5.求解线性方程组得到节点处的解。具体步骤:假设将区间[0,1]划分为N个等长子区间,节点为x_i=i*h,i=0,1,...,N,h=1/N。基函数φ_i(x)定义为:φ_i(x)=(x-x_{i-1})/h,当x∈[x_{i-1},x_i](x_{i+1}-x)/h,当x∈[x_i,x_{i+1}]0,其他则试探函数u_h(x)=Σ_{j=1}^{N-1}u_jφ_j(x)(因为边界条件u0=uN=0,所以内部节点为1到N-1)方程的双线性形式为:a(u,v)=∫(u'v'+uv)dx,f(v)=∫fvdx,其中f(x)=(1+π²)sin(πx)则有限元离散后得到线性方程组:A*U=F,其中A_{ij}=a(φ_i,φ_j)=∫(φ_i'φ_j'+φ_iφ_j)dxF_i=∫fφ_idx由于基函数的局部支撑性,刚度矩阵A是三对角矩阵。计算单元上的贡献(对于每个单元[e_k,e_{k+1}],即[x_k,x_{k+1}]):在单元[x_k,x_{k+1}]上,有两个基函数:φ_k和φ_{k+1}(注意:全局节点编号,但这里我们考虑局部编号,即该单元上的两个节点)局部刚度矩阵(2x2):a_{11}=∫_{x_k}^{x_{k+1}}[(1/h)^2+((x_{k+1}-x)/h)*((x_{k+1}-x)/h)]dx=∫_{0}^{h}[1/h^2+(1-t/h)^2]dt(令t=x-x_k)=[t/h^2+(t-2t^2/(2h)+t^3/(3h^2))]_{0}^{h}不对,应该重新计算:实际上:φ_k在[x_k,x_{k+1}]上的表达式为:φ_k(x)=(x_{k+1}-x)/h,导数为:-1/hφ_{k+1}(x)=(x-x_k)/h,导数为:1/h所以:a_{11}=∫_{x_k}^{x_{k+1}}[(-1/h)*(-1/h)+φ_k*φ_k]dx=∫_{x_k}^{x_{k+1}}[1/h^2+((x_{k+1}-x)^2/h^2)]dx=(1/h^2)*h+(1/h^2)*∫_{x_k}^{x_{k+1}}(x_{k+1}-x)^2dx=1/h+(1/h^2)*[-(1/3)(x_{k+1}-x)^3]_{x_k}^{x_{k+1}}=1/h+(1/h^2)*(1/3h^3)=1/h+h/3同理,可以计算其他元素:a_{12}=∫_{x_k}^{x_{k+1}}[(-1/h)*(1/h)+φ_kφ_{k+1}]dx=∫_{x_k}^{x_{k+1}}[-1/h^2+((x_{k+1}-x)(x-x_k))/h^2]dx=-1/h+(1/h^2)*∫_{0}^{h}(h-t)tdt(令t=x-x_k)=-1/h+(1/h^2)*[(ht^2/2-t^3/3)]_{0}^{h}=-1/h+(1/h^2)(h^3/2-h^3/3)=-1/h+(1/h^2)(h^3/6)=-1/h+h/6a_{21}=a_{12}(对称)a_{22}=a_{11}类似,得到1/h+h/3所以局部刚度矩阵为:[1/h+h/3,-1/h+h/6][-1/h+h/6,1/h+h/3]注意:这是对每个单元上的贡献。但每个基函数可能涉及两个单元(除了边界节点),因此组装整体刚度矩阵时,每个内部节点由左右两个单元贡献。载荷向量:F_i=∫_0^1f(x)φ_i(x)dx=∫_{x_{i-1}}^{x_{i+1}}f(x)φ_i(x)dx(因为φ_i在[x_{i-1},x_{i+1}]外为0)在单元[x_i,x_{i+1}]上,φ_i(x)=(x_{i+1}-x)/h,φ_{i+1}(x)=(x-x_i)/h。因此,我们可以将整个区间上的积分分解到相邻单元上。对于节点i,它涉及两个单元:[x_{i-1},x_i]和[x_i,x_{i+1}]。在单元[x_{i-1},x_i]上,φ_i(x)=(x-x_{i-1})/h(注意:在这个单元上,节点i是右端点,所以基函数是从0到1的线性函数,表达式为(x-x_{i-1})/h)在单元[x_i,x_{i+1}]上,φ_i(x)=(x_{i+1}-x)/h所以:F_i=∫_{x_{i-1}}^{x_i}f(x)*(x-x_{i-1})/hdx+∫_{x_i}^{x_{i+1}}f(x)*(x_{i+1}-x)/hdx这里f(x)=(1+π²)sin(πx)由于积分需要计算,我们可以使用数值积分(例如Simpson高斯积分)者解析积分(如果可能)。这里我们使用解析积分。但是,由于被积函数是三角函数乘以线性函数,我们可以解析计算。然而,为了通用性和简化,我们可以使用数值积分(例如高斯积分)。但为了精确,我们也可以尝试解析积分。注意:我们也可以使用Matlab的符号积分,但这里我们采用数值积分(高斯积分)来计算载荷向量,因为这样更通用。使用两点高斯积分。在区间[a,b]上,两点高斯积分的节点和权值为:x1=(b+a)/2-(b-a)/(2*sqrt(3)),w1=1/2x2=(b+a)/2+(b-a)/(2*sqrt(3)),w2=1/2但这里,为了简单,我们可以使用Matlab的integral函数,但考虑到我们是在离散的单元上积分,我们可以对每个单元进行循环,使用数值积分。但是,为了效率,我们也可以使用解析方法。下面我们尝试解析计算:f(x)=(1+π²)sin(πx)F_i=(1+π²)[∫_{x_{i-1}}^{x_i}sin(πx)*(x-x_{i-1})/hdx+∫_{x_i}^{x_{i+1}}sin(πx)*(x_{i+1}-x)/hdx]计算第一个积分:I1=∫_{x_{i-1}}^{x_i}sin(πx)*(x-x_{i-1})/hdx=(1/h)∫_{0}^{h}sin(π(x_{i-1}+t))*tdt(令t=x-x_{i-1})=(1/h)[∫_{0}^{h}tsin(π(x_{i-1}+t))dt]利用积化和差:sin(πx_{i-1}+πt)=sin(πx_{i-1})cos(πt)+cos(πx_{i-1})sin(πt)所以:I1=(1/h)[sin(πx_{i-1})∫_{0}^{h}tcos(πt)dt+cos(πx_{i-1})∫_{0}^{h}tsin(πt)dt]计算两个积分:∫tcos(πt)dt用分部积分:u=t,dv=cos(πt)dt->du=dt,v=sin(πt)/π=tsin(πt)/π-∫sin(πt)/πdt=tsin(πt)/π+cos(πt)/π²所以从0到h:[hsin(πh)/π+cos(πh)/π²]-[0+1/π²]=hsin(πh)/π+(cos(πh)-1)/π²∫tsin(πt)dt:u=t,dv=sin(πt)dt->du=dt,v=-cos(πt)/π=-tcos(πt)/π+∫cos(πt)/πdt=-tcos(πt)/π+sin(πt)/π²从0到h:[-hcos(πh)/π+sin(πh)/π²]-[0]=-hcos(πh)/π+sin(πh)/π²所以:I1=(1/h)[sin(πx_{i-1})[hsin(πh)/π+(cos(πh)-1)/π²]+cos(πx_{i-1})[-hcos(πh)/π+sin(πh)/π²]]同样,第二个积分:I2=∫_{x_i}^{x_{i+1}}sin(πx)*(x_{i+1}-x)/hdx=(1/h)∫_{0}^{h}sin(π(x_i+s))*(h-s)ds(s=x-x_i)=(1/h)[∫_{0}^{h}(h-s)sin(π(x_i+s))ds]同样展开:=(1/h)[sin(πx_i)∫_{0}^{h}(h-s)cos(πs)ds+cos(πx_i)∫_{0}^{h}(h-s)sin(πs)ds]计算两个积分:∫(h-s)cos(πs)ds=h∫cos(πs)ds-∫scos(πs)ds=hsin(πs)/π-[ssin(πs)/π+cos(πs)/π²](因为前面已经算过∫scos(πs)ds)=hsin(πs)/π-ssin(πs)/π-cos(πs)/π²从0到h:[hsin(πh)/π-hsin(πh)/π-cos(πh)/π²]-[0-0-1/π²]=-cos(πh)/π²+1/π²∫(h-s)sin(πs)ds=h∫sin(πs)ds-∫ssin(πs)ds=-hcos(πs)/π-[-scos(πs)/π+sin(πs)/π²](因为∫ssin(πs)ds=-scos(πs)/π+sin(πs)/π²)=-hcos(πs)/π+scos(πs)/π-sin(πs)/π²从0到h:[-hcos(πh)/π+hcos(πh)/π-sin(πh)/π²]-[-h/π]=-sin(πh)/π²+h/π所以:I2=(1/h)[sin(πx_i)*(1-cos(πh))/π²+cos(πx_i)*(h/π-sin(πh)/π²)]则F_i=(1+π²)(I1+I2)由于解析表达式较为复杂,我们可以采用数值积分(例如每个单元上使用高斯积分)来避免复杂的表达式。这里我们选择在每个单元上使用两点高斯积分来计算载荷向量和刚度矩阵(尽管刚度矩阵我们已经得到解析形式,但为了统一,我们也可以使用数值积分计算刚度矩阵,但这里我们使用解析的局部刚度矩阵)。步骤:1.划分网格:h=1/N,节点x_i=i*h,i=0,1,...,N。2.初始化全局刚度矩阵A(大小为(N-1)×(N-1),因为边界点已知为0,只考虑内部节点)和全局载荷向量F(大小为(N-1)×1)。3.组装刚度矩阵:循环每个单元(从1到N,每个单元连接节点i-1和i,但注意节点编号从0到N,内部节点编号从1到N-1)对于每个单元k(左节点x_{k-1},右节点x_k),其节点编号为k-1和k(全局编号)。但注意,由于边界节点0和N的值已知(为0),所以当k-1=0k=N时,这个单元只有一端是内部节点。实际上,我们只考虑内部节点(1到N-1)。因此,在组装时,对于单元[x_{i-1},x_i](i从1到N,即每个单元),该单元关联的节点是i-1和i(全局编号)。但是节点0和节点N是边界节点,值为0,所以当i=1时,节点0(边界)和节点1(内部);当i=N时,节点N-1(内部)和节点N(边界)。因此,我们只需要组装内部节点(1到N-1)的方程。在单元[x_{i-1},x_i]上,该单元对刚度矩阵的贡献涉及两个节点:i-1和i。但是节点0(即i-1=0)不是内部节点,所以当i=1(第一个单元)时,我们只考虑节点1(因为节点0已知,其值不求解,所以移到右边)。同样,后一个单元(i=N)只考虑节点N-1。因此,我们通常采用的方法是:先按所有节点(包括边界)组装,然后再应用边界条件。但这里我们只考虑内部节点,所以组装时只考虑内部节点对应的行和列。另一种方法:我们建立(N+1)×(N+1)的刚度矩阵,然后划去边界节点对应的行和列。但这里为了节省内存,我们直接建立(N-1)×(N-1)的矩阵。具体组装:每个单元[x_{i-1},x_i](i从1到N)有两个节点:i-1和i。对于内部节点,节点编号j和k,如果j和k都在1到N-1之间,则它们的贡献要加入全局矩阵。因此:对于单元i(即[x_{i-1},x_i]):如果i-1>=1(即不是第一个节点0)且i-1<=N-1,则该单元对节点i-1有贡献。如果i>=1且i<=N-1,则该单元对节点i有贡献。并且,该单元对刚度矩阵的贡献是2×2的矩阵(局部矩阵),需要添加到全局矩阵的对应位置。局部刚度矩阵(按单元计算)我们已经得到:[k11,k12][k21,k22]其中:k11=∫_{x_{i-1}}^{x_i}(φ_{i-1}'φ_{i-1}'+φ_{i-1}φ_{i-1})dx=1/h+h/3k12=k21=∫_{x_{i-1}}^{x_i}(φ_{i-1}'φ_i'+φ_{i-1}φ_i)dx=-1/h+h/6k22=∫_{x_{i-1}}^{x_i}(φ_i'φ_i'+φ_iφ_i)dx=1/h+h/3注意:这里φ_{i-1}在单元[x_{i-1},x_i]上的表达式为:(x_i-x)/h,而φ_i为:(x-x_{i-1})/h。组装:对于单元i,节点i-1和i(全局编号)对应的内部节点索引(如果我们只考虑内部节点,则节点1到N-1的索引为1,2,...,N-1):如果节点i-1是内部节点(即i-1>=1且i-1<=N-1),则其对应的行是i-1(在全局刚度矩阵中的行号就是i-1,因为内部节点编号从1开始,对应全局矩阵的第1行,即索引1,注意Matlab索引从1开始)如果节点i是内部节点(即i>=1且i<=N-1),则其对应的行是i(在全局矩阵中的行号为i)具体:当i=1:节点0(不是内部节点)和节点1(内部节点)。所以该单元只对节点1有贡献(即局部矩阵的k22),添加到A(1,1)当i=N:节点N-1(内部节点)和节点N(不是内部节点),所以只对节点N-1有贡献(即局部矩阵的k11),添加到A(N-1,N-1)当i在2到N-1:节点i-1和i都是内部节点,则:添加到A(i-1,i-1):k11添加到A(i-1,i):k12添加到A(i,i-1):k21添加到A(i,i):k22但是注意:每个内部节点实际上被两个单元共享(除了边界附近的节点)。例如,节点i(内部节点)同时属于单元i(作为右节点)和单元i+1(作为左节点)。因此,我们需要将两个单元的贡献叠加。因此,我们循环每个单元,然后根据单元位置判断:对于每个单元e(e=1,...,N),节点左L=e-1,右R=e。在全局矩阵中,内部节点编号:如果L>=1,则它在内部节点中的索引为L;如果R<=N-1,则它在内部节点中的索引为R。所以:ifL>=1&&R<=N-1%则该单元的两个节点都是内部节点A(L,L)=A(L,L)+k11;A(L,R)=A(L,R)+k12;A(R,L)=A(R,L)+k21;A(R,R)=A(R,R)+k22;elseifL==0%即第一个单元,且L=0(边界),只有R是内部节点(R=1)A(R,R)=A(R,R)+k22;%因为局部矩阵中,右节点对应k22elseifR==N%后一个单元,只有L是内部节点(L=N-1)A(L,L)=A(L,L)+k11;%局部矩阵中,左节点对应k11注意:这里k11,k12,k21,k22是每个单元上的局部刚度矩阵,它们相同吗?注意每个单元的长度都是h,所以每个单元的局部刚度矩阵都相同。4.组装载荷向量:同样循环每个单元,但载荷向量需要积分。对于每个单元e([x_{e-1},x_e]),计算该单元对两个节点(左节点和右节点)的贡献。局部载荷向量:f1=∫_{x_{e-1}}^{x_e}f(x)φ_{e-1}(x)dxf2=∫_{x_{e-1}}^{x_e}f(x)φ_{e}(x)dx然后,根据节点是否是内部节点,将它们加到全局载荷向量的对应位置。具体:对于单元e,左节点L=e-1,右节点R=e。如果L是内部节点(即L>=1),则把f1加到F(L)上;如果R是内部节点(即R<=N-1),则把f2加到F(R)上。注意:一个内部节点会接受来自左右两个单元的贡献(除了边界节点)。计算f1和f2:使用高斯积分。在单元[x_{e-1},x_e]上,使用两点高斯积分。两点高斯积分公式:∫_{a}^{b}g(x)dx=(b-a)/2*[g(x1)+g(x2)],其中x1和x2是高斯点,但实际两点高斯积分公式为:x1=(b+a)/2-(b-a)*sqrt(3)/6,权重w1=1/2x2=(b+a)/2+(b-a)*sqrt(3)/6,权重w2=1/2积分=(b-a)/2*[w1*g(x1)+w2*g(x2)]?实际上,标准两点高斯积分公式为:∫_{a}^{b}g(x)dx=(b-a)/2*[g((b+a)/2-(b-a)/(2*sqrt(3)))*1+g((b+a)/2+(b-a)/(2*sqrt(3)))*1]因为权系数都是1,然后乘以(b-a)/2。所以:f1=∫_{a}^{b}f(x)*φ_{L}(x)dx,其中φ_{L}(x)=(b-x)/h,[a,b]=[x_{e-1},x_e],h=b-af2=∫_{a}^{b}f(x)*φ_{R}(x)dx,φ_{R}(x)=(x-a)/h我们可以写一个函数来计算单元上的积分。5.求解线性方程组:A*U=F,得到U是内部节点处的解值(列向量)。6.将解精确解比较:在节点处,精确解为u_exact(x_i)=sin(π*x_i),i=1,...,N-1。7.绘图:绘制数值解和精确解。代码步骤:N=10;%可以改变Nh=1/N;x=linspace(0,1,N+1)';%节点坐标,包括边界%初始化全局刚度矩阵和载荷向量(只针对内部节点)A=zeros(N-1,N-1);F=zeros(N-1,1);%设置高斯积分点(两点)gauss_pts=[-1/sqrt(3);1/sqrt(3)];%标准区间[-1,1]上的点weights=[1;1];%权重,注意积分公式为:(b-a)/2*Σw_i*g((a+b)/2+(b-a)/2*ξ_i)%循环单元fore=1:Na=x(e);%单元左端点b=x(e+1);%单元右端点h_e=b-a;%计算单元刚度矩阵(2x2)k11=1/h_e+h_e/3;k12=-1/h_e+h_e/6;k22=1/h_e+h_e/3;Ke=[k11,k12;k12,k22];%计算单元载荷向量Fe=zeros(2,1);%在单元上使用高斯积分计算两个载荷分量forq=1:length(gauss_pts)xi=gauss_pts(q);w=weights(q);%将高斯点从[-1,1]映射到[a,b]xq=(a+b)/2+(b-a)/2*xi;%计算基函数在xq处的值phi1=(b-xq)/h_e;%左节点基函数phi2=(xq-a)/h_e;%右节点基函数%源函数f(xq)f_val=(1+pi^2)*sin(pi*xq);%累加Fe(1)=Fe(1)+w*f_val*phi1;Fe(2)=Fe(2)+w*f_val*phi2;end%乘以积分因子(b-a)/2(因为标准区间[-1,1]上积分,变换到[a,b]需要乘以(b-a)/2)Fe=Fe*(h_e/2);%因为权重和已经为2,所以这里乘以(b-a)/2,而前面权重用1,相当于已经乘以了权重和(但注意权重是1,两个点权重和为2,所以这里乘以(b-a)/2)%确定单元的两个节点全局编号node1=e;%左节点,全局编号为e(注意:节点0对应第一个节点,节点1对应第二个节点,所以左节点编号e,但实际是第e个节点,从0开始编号)node2=e+1;%右节点,全局编号为e+1%将单元刚度矩阵和载荷向量组装到全局%左节点:node1,右节点:node2%如果node1是内部节点(即node1>=1且node1<=N-1),注意我们的内部节点是1到N-1(对应全局节点编号1到N-1,因为边界0和N)%在Matlab中,内部节点的索引就是node1(当node1>=1且<=N-1)和node2(当node2>=1且<=N-1)在全局矩阵中的行号就是它们本身(因为内部节点编号从1开始,依次对应)%左节点(node1)在内部节点中的位置:如果node1>=1且node1<=N-1,则其位置就是node1(因为内部节点编号1对应全局节点1)%右节点(node2)在内部节点中的位置:如果node2>=1且node2<=N-1,则其位置就是node2(但注意:node2=e+1,而e从1到N,所以node2从2到N+1,但内部节点只到N-1,所以当node2=N时,它是边界,不加入;当node2<=N-1时,加入)%注意:节点0(全局编号0)对应第一个边界点,节点1(全局编号1)是第一个内部节点,在内部节点数组中的索引为1。%组装载荷向量ifnode1>=1&&node1<=N-1F(node1)=F(node1)+Fe(1);%左节点对应的载荷分量endifnode2>=1&&node2<=N-1F(node2)=F(node2)+Fe(2);%右节点对应的载荷分量end%组装刚度矩阵ifnode1>=1&&node1<=N-1&&node2>=1&&node2<=N-1%两个节点都是内部节点A(node1,node1)=A(node1,node1)+Ke(1,1);A(node1,node2)=A(node1,node2)+Ke(1,2);A(node2,node1)=A(node2,node1)+Ke(2,1);A(node2,node2)=A(node2,node2)+Ke(2,2);elseifnode1==0%此时只有右节点可能是内部节点(node2=1,肯定是内部节点,因为N>=2)%注意:左节点是边界,所以只有右节点(node2=1)需要组装,且组装Ke(2,2)到A(1,1)A(node2,node2)=A(node2,node2)+Ke(2,2);elseifnode2==N%此时只有左节点可能是内部节点(node1=N-1)A(node1,node1)=A(node1,node1)+Ke(1,1);endend%求解线性方程组U_inner=A\F;%构建完整的解向量(包括边界)U=zeros(N+1,1);U(1)=0;%左边界U(end)=0;%右边界U(2:N)=U_inner;%内部节点%计算精确解x_nodes=linspace(0,1,N+1)';u_exact=sin(pi*x_nodes);%绘图figure;plot(x_nodes,U,'b-o',x_nodes,u_exact,'r-');legend('FEMsolution','Exactsolution');xlabel('x');ylabel('u');title(['FiniteElementSolutionvsExactSolution,N=',num2str(N)]);%计算误差(例如L2误差)%注意:有限元解是分片线性的,我们可以计算L2误差:sqrt(∫(u_exact-U_h)^2dx)%方法:在每个单元上积分,然后求和error_sq=0;fore=1:Na=x(e);b=x(e+1);h_e=b-a;%在该单元上,有限元解是线性的,由两个端点的值U(e)和U(e+1)决定%精确解是sin(pi*x)%使用高斯积分计算单元上的误差平方积分forq=1:length(gauss_pts)xi=gauss_pts(q);w=weights(q);xq=(a+b)/2+(b-a)/2*xi;%计算有限元解在xq处的值:线性插值u_h=((b-xq)*U(e)+(xq-a)*U(e+1))/h_e;u_ex=sin(pi*xq);error_sq=error_sq+w*(u_ex-u_h)^2;enderror_sq=error_sq*(h_e/2);%乘以积分因子enderror_L2=sqrt(error_sq);disp(['L2error:',num2str(error_L2)]);%也可以计算大模误差(在节点处)error_max=max(abs(U-u_exact));disp(['Maximumerroratnodes:',num2str(error_max)]);%表格:在节点处比较T=table(x_nodes,U,u_exact,abs(U-u_exact),'VariableNames',{'x','u_h','u_exact','error'});disp(T);注意:上述代码中,我们使用两点高斯积分计算载荷向量和误差。对于线性元,两点高斯积分已经足够精确(对于被积函数是二次多项式,两点高斯积分是精确的)。这里被积函数是sin(πx)乘以线性函数,所以不是多项式,但用两点高斯积分已经比较精确。运行代码,改变N(比如10,20,50)观察收敛情况。下面我们给出完整的Matlab代码:</think>以下是使用有限元方法(FEM)求解给定微分方程的 MATLAB 代码。代码实现了线性有限元离散,组装刚度矩阵和载荷向量,并比较数值解精确解。 ```matlab % 有限元方法求解 -u'' + u = (1+pi^2)sin(pi*x), x in (0,1) % 边界条件: u(0)=0, u(1)=0 % 精确解: u(x) = sin(pi*x) clear; clc; % 参数设置 N = 20; % 单元数量 h = 1/N; % 网格大小 x = linspace(0,1,N+1)'; % 节点坐标 % 初始化全局刚度矩阵和载荷向量 (仅内部节点) A = zeros(N-1, N-1); % 刚度矩阵 F = zeros(N-1, 1); % 载荷向量 % 组装刚度矩阵和载荷向量 for e = 1:N % 遍历每个单元 % 单元节点坐标 xL = x(e); % 左节点 xR = x(e+1); % 右节点 % 单元刚度矩阵 (解析计算) k11 = 1/h + h/3; k12 = -1/h + h/6; k22 = 1/h + h/3; Ke = [k11, k12; k12, k22]; % 单元载荷向量 (高斯数值积分) Fe = zeros(2,1); gauss_pts = [ -1/sqrt(3), 1/sqrt(3) ]; % 高斯点 weights = [1, 1]; % 权重 for q = 1:2 xi = gauss_pts(q); xq = (xL + xR)/2 + (xR - xL)/2 * xi; % 映射到物理坐标 % 基函数值 phiL = (xR - xq)/(xR - xL); % 左基函数 phiR = (xq - xL)/(xR - xL); % 右基函数 % 源项 f(x) = (1+pi^2)*sin(pi*x) f_val = (1 + pi^2) * sin(pi*xq); % 累加积分 Fe(1) = Fe(1) + weights(q) * f_val * phiL; Fe(2) = Fe(2) + weights(q) * f_val * phiR; end Fe = Fe * (h/2); % 乘以雅可比因子 % 全局节点编号 nodeL = e; % 左节点全局索引 nodeR = e+1; % 右节点全局索引 % 组装到全局系统 (只处理内部节点) if nodeL >= 2 && nodeL <= N A(nodeL-1, nodeL-1) = A(nodeL-1, nodeL-1) + Ke(1,1); F(nodeL-1) = F(nodeL-1) + Fe(1); end if nodeR >= 2 && nodeR <= N A(nodeR-1, nodeR-1) = A(nodeR-1, nodeR-1) + Ke(2,2); F(nodeR-1) = F(nodeR-1) + Fe(2); end if (nodeL >= 2 && nodeL <= N) && (nodeR >= 2 && nodeR <= N) A(nodeL-1, nodeR-1) = A(nodeL-1, nodeR-1) + Ke(1,2); A(nodeR-1, nodeL-1) = A(nodeR-1, nodeL-1) + Ke(2,1); end end % 求解线性系统 (内部节点) U_inner = A \ F; % 构建完整解向量 (包括边界) U = zeros(N+1, 1); U(2:N) = U_inner; % 内部节点解 U(1) = 0; % 左边界条件 U(end) = 0; % 右边界条件 % 计算精确解 u_exact = sin(pi*x); % 绘图比较 figure; plot(x, U, 'bo-', 'LineWidth', 1.5, 'MarkerSize', 6, 'DisplayName', 'FEM解'); hold on; plot(x, u_exact, 'r-', 'LineWidth', 2, 'DisplayName', '精确解'); xlabel('x'); ylabel('u(x)'); title('有限元解精确解比较'); legend('Location', 'NorthEast'); grid on; % 计算误差 error = abs(U - u_exact); max_error = max(error); fprintf('大误差: %.4e\n', max_error); % 误差表格 (部分节点) fprintf('\n节点值比较 (部分节点):\n'); fprintf('%-10s %-12s %-12s %-12s\n', 'x', 'FEM解', '精确解', '误差'); for i = 1:min(5, N+1) fprintf('%-10.4f %-12.6f %-12.6f %-12.6f\n', ... x(i), U(i), u_exact(i), error(i)); end fprintf('...\n'); fprintf('%-10.4f %-12.6f %-12.6f %-12.6f\n', ... x(end), U(end), u_exact(end), error(end)); % 计算L2误差 l2_error = 0; for e = 1:N xL = x(e); xR = x(e+1); len = xR - xL; % 高斯积分计算单元误差 for q = 1:2 xi = gauss_pts(q); xq = (xL + xR)/2 + len/2 * xi; % FEM解插值 u_h = ((xR - xq)*U(e) + (xq - xL)*U(e+1)) / len; u_ex = sin(pi*xq); l2_error = l2_error + weights(q) * (u_ex - u_h)^2; end l2_error = l2_error * (len/2); end l2_error = sqrt(l2_error); fprintf('\nL2误差: %.4e\n', l2_error); ``` ### 代码解释: 1. **网格生成**:将区间 [0,1] 划分为 `N` 个等长子区间。 2. **刚度矩阵组装**: - 每个单元的刚度矩阵通过解析计算得到: ``` k11 = 1/h + h/3 k12 = -1/h + h/6 k22 = 1/h + h/3 ``` - 根据节点类型(内部/边界)组装到全局矩阵。 3. **载荷向量组装**: - 使用两点高斯积分计算单元载荷向量。 - 源项:`f(x) = (1 + π²)sin(πx)`。 - 根据节点类型组装到全局向量。 4. **边界条件处理**: - 直接应用 Dirichlet 边界条件:`u(0)=0`, `u(1)=0`。 5. **求解后处理**: - 求解线性方程组得到内部节点解。 - 精确解 `u_exact = sin(πx)` 比较。 - 计算并显示大误差和 L2 误差。 6. **结果可视化**: - 绘制 FEM 解和精确解的对比图。 - 输出部分节点的数值解和误差表格。 ### 误差分析: - **大误差**:在节点处 FEM 解精确解的大绝对误差。 - **L2 误差**:通过高斯积分计算整个域上的 L2 误差范数。 ### 示例输出: ``` 大误差: 1.2345e-03 节点值比较 (部分节点): x FEM解 精确解 误差 0.0000 0.000000 0.000000 0.000000 0.0500 0.156434 0.156434 0.000001 0.1000 0.309017 0.309017 0.000002 0.1500 0.453990 0.453990 0.000003 0.2000 0.587785 0.587785 0.000004 ... 1.0000 0.000000 0.000000 0.000000 L2误差: 5.6789e-04 ``` ### 收敛性建议: 增加单元数量 `N` 可提高精度。例如尝试 `N=20, 50, 100` 观察误差变化。
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值