报告作业报告大作业报告

武汉大学数值实习MATLAB作业期末报告

时间：2024.4.20

实  验  报  告

实验课程：MATLAB数值分析

学生姓名：        

学    号： 20103010

学    院：

专    业：     

             指导老师：         

20##年12月28日

一、问题描述

考虑如下Dirichlet问题

  其中为正方形区域的边界，类似于模型问题，我们得到差分方程



得到系数矩阵

A=其中B=-I，I为n-1阶单位矩阵S’为对角元为1+次对角元为的n-1阶的对称三对角矩阵。对=sin（xy），=，n=20.用超松弛以及共轭梯度求解差分方程的解。

二、Matlab代码及其结果

Gauss-Seidel方法：

n=20;

h=1/n;

a=h*h/4;

m=n*n;

S=zeros(n-1,n-1);

for i=1:n-1

    S(i,i)=1+h*h/4;

end

for i=2:n-1

    S(i,i-1)=-1/4;

end

for i=1:n-2

    S(i,i+1)=-1/4;

end

B=zeros(n-1,n-1);

for i=1:n-1

    B(i,i)=-1/4;

end

A=zeros((n-1)*(n-1),(n-1)*(n-1));

for i=1:n-1

    A((i-1)*(n-1)+1:i*(n-1),(i-1)*(n-1)+1:i*(n-1))=S;

end

for i=1:n-2

    A(i*(n-1)+1:(i+1)*(n-1),(i-1)*(n-1)+1:i*(n-1))=B;

end

for i=1:n-2

    A((i-1)*(n-1)+1:i*(n-1),i*(n-1)+1:(i+1)*(n-1))=B;

End

b=[1:(n-1)*(n-1)]'*0;

j=1;                     

b(1)=a*sin(1/m)+2*a;

b(n-1)=a*sin((n-1)/m)+a*((n-1)*(n-1))+a*(1+n*n);

for i=2:n-2

    b(i)=a*sin(i*j/m)+a*i*i;

end

for j=2:n-2

    b((j-1)*(n-1)+1)=a*sin(1*j/m)+a*j*j;

    b(j*(n-1))=a*sin((n-1)*j/m)+a*(j*j+n*n);

    for i=(j-1)*(n-1)+2:j*(n-1)-1

         b(i)=a*sin((i-(j-1)*(n-1))*j/m);

    end

end

j=n-1;

b((n-2)*(n-1)+1)=a*sin(1*j/m)+(n-1)*(n-1)*a+(1+n*n)*a;

for i=(n-2)*(n-1)+2:(n-1)*(n-1)-1

    b(i)=a*sin((i-(n-2)*(n-1))*j/m)+1/4+((i-(n-2)*(n-1))*(n-1))*a;

end

b((n-1)*(n-1))=a*sin((n-1)*(n-1)/m)+(n*n+(n-1)*(n-1))/(m*2);

y=A\b

n=length (b);

Y=zeros (n,1);

Y0=Y;

Y1=Y;

w=1.7

U=zeros (n,n);

for i=1:n

    for j=i+1:n

        U (i,j)=A (i,j);

    end

end

DL=A-U;

for II=1:100000

    B=-U*Y0+b;

    Y1 (1)=B (1)/DL(1,1);

    for i=2:n

        Y1 (i)=(B (i)-DL (i,1:i-1)*Y1 (1:i-1))/DL(i,i);

    end

    Y=w*Y1+(1-w)*Y0;

    Y0=Y1;

   

    r=A*Y-b;

    err=norm (r);

    if err<eps

        break;

    end

end

[II,err]

Y

reshape(Y,19,19)

SOR方法：

n=20;

h=1/n;

a=h*h/4;

m=n*n;

S=zeros(n-1,n-1);

for i=1:n-1

    S(i,i)=1+h*h/4;

end

for i=2:n-1

    S(i,i-1)=-1/4;

end

for i=1:n-2

    S(i,i+1)=-1/4;

end

B=zeros(n-1,n-1);

for i=1:n-1

    B(i,i)=-1/4;

end

A=zeros((n-1)*(n-1),(n-1)*(n-1));

for i=1:n-1

    A((i-1)*(n-1)+1:i*(n-1),(i-1)*(n-1)+1:i*(n-1))=S;

end

for i=1:n-2

    A(i*(n-1)+1:(i+1)*(n-1),(i-1)*(n-1)+1:i*(n-1))=B;

end

for i=1:n-2

    A((i-1)*(n-1)+1:i*(n-1),i*(n-1)+1:(i+1)*(n-1))=B;

end

b=[1:(n-1)*(n-1)]'*0;

j=1;                     

b(1)=a*sin(1/m)+2*a;

b(n-1)=a*sin((n-1)/m)+a*((n-1)*(n-1))+a*(1+n*n);

for i=2:n-2

    b(i)=a*sin(i*j/m)+a*i*i;

end

for j=2:n-2

    b((j-1)*(n-1)+1)=a*sin(1*j/m)+a*j*j;

    b(j*(n-1))=a*sin((n-1)*j/m)+a*(j*j+n*n);

    for i=(j-1)*(n-1)+2:j*(n-1)-1

         b(i)=a*sin((i-(j-1)*(n-1))*j/m);

    end

end

j=n-1;

b((n-2)*(n-1)+1)=a*sin(1*j/m)+(n-1)*(n-1)*a+(1+n*n)*a;

for i=(n-2)*(n-1)+2:(n-1)*(n-1)-1

    b(i)=a*sin((i-(n-2)*(n-1))*j/m)+1/4+((i-(n-2)*(n-1))*(n-1))*a;

end

b((n-1)*(n-1))=a*sin((n-1)*(n-1)/m)+(n*n+(n-1)*(n-1))/(m*2);

y=A\b；

Y=b*0

p=Y;

r=Y;

p1=Y;

r1=Y;

r=b-A*Y;

k=0;

normr=r'*r;

while sqrt (normr)>eps

    k=k+1

    if k==1

        p=r

    else

        beta=normr2/normr;

        p=r+beta*p;

        normr=normr2;

    end

    alpha=(normr)/ (p'*A* p);

    Y=Y+alpha*p;

    r=r-alpha*A*p;

    normr2=r'*r;

    k;

    normr2;

end

sqrt_normr=sqrt (normr2)

k

r

Y 

reshape(Y,19,19)


结果图片： 

    

三、C语言代码及其结果

Gauss-Seidel方法：

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

int main()

{

      int i,j,k;

      int m,n,p;

      double b[400]={0},A[400][400]={0},U[400]={0};

      double P[400]={0},r[400]={0},B[400]={0};

      double a;

      double h;

      double eps;

      double normr=0.0,normr2,alpha,beta;

      double sum,sumr,sumA,sqrt_normr;

      printf("n= \n");

      scanf("%d",&n);

      printf("eps= \n");

      scanf("%lf",&eps);

      printf("%lf\n",eps);

      h=1.0/n;

      m=n*n;

      p=n-1;

      printf("%lf\n",h);

      a=h*h/4;

      printf("%d,%d,%lf,%lf\n",m,p,h,a);

      j=1;

      b[1]=a*sin(1.0/m)+2*a;

      b[n-1]=a*sin(1.0*p/m)+a*(p*p)+a*(1+n*n);

      for(i=2;i<=n-2;i++)

          b[i]=a*sin(1.0*i*j/m)+a*i*i;

      for(j=2;j<=n-2;j++)

      {

          b[(j-1)*p+1]=a*sin(1.0*j/m)+a*j*j;

          b[j*p]=a*sin(1.0*p*j/m)+a*(j*j+n*n);

          for(i=(j-1)*p+2;i<=j*p-1;i++)

              b[i]=a*sin((i-(j-1)*p)*j*1.0/m);

      }

          /*for(i=17*p+2;i<=18*p;i++)

         printf("%0.10lf\n",b[i]);*/

           j=n-1;

      b[(p-1)*p+1]=a*sin(1.0*j/m)+p*p*a+a*(1+n*n);

           for(i=(p-1)*p+2;i<=p*p-1;i++)

          b[i]=a*sin(1.0*(i-(p-1)*p)*j/m)+1.0/4+a*(i-(p-1)*p)*(i-(p-1)*p);

      b[p*p]=a*sin(1.0*p*p/m)+(n*n+p*p)*2*a;

        A[1][1]=1+a;

      A[1][2]=-0.25;

      A[p*p][p*p-1]=-0.25;

      A[p*p][p*p]=1+a;

      for(i=2;i<=p*p-1;i++)

      {   A[i][i-1]=-0.25;

          A[i][i]=1+a;

          A[i][i+1]=-0.25;

      }

      for(i=n;i<=p*p;i++)

      {    A[i-n+1][i]=-0.25;

           A[i][i-n+1]=-0.25;

      }

      for(j=2;j<=n-1;j++)

      {   A[(j-1)*p][(j-1)*p+1]=0.0;

          A[(j-1)*p+1][(j-1)*p]=0.0;

      }

                for(i=1;i<=p*p;i++)

        {    sum=0.0;

             for(j=1;j<=p*p;j++)

                 sum=sum+A[i][j]*U[j];

             r[i]=b[i]-sum;

        }

                   for(i=1;i<=p*p;i++)

            normr=normr+r[i]*r[i];

             k=0;

        do

        {  

            k=k+1;

            if(k==1)

                {

                    for(i=1;i<=p*p;i++)

                        P[i]=r[i];                    

                }

            else

                {

                    beta=normr2/normr;

                    for(i=1;i<=p*p;i++)

                        P[i]=r[i]+beta*P[i];

                    normr=normr2;              

                 }

            sumr=0.0;     //r的范数//

            for(i=1;i<=p*p;i++)

                sumr=sumr+r[i]*r[i];

                     for(i=1;i<=p*p;i++)      //求A*P//

              {   

                   B[i]=0.0;

                   for(j=1;j<=p*p;j++)

                      B[i]=B[i]+A[i][j]*P[j];

              }

                     sumA=0.0;      //A关于P的共轭乘积//

            for(i=1;i<=p*p;i++)

               sumA=sumA+P[i]*B[i];

                      alpha=sumr/sumA;     //求alpha//

                      for(i=1;i<=p*p;i++)

                U[i]=U[i]+alpha*P[i];

            for(i=1;i<=p*p;i++)

                r[i]=r[i]-alpha*B[i];

                    normr2=0.0;

            for(i=1;i<=p*p;i++)

                normr2=normr2+r[i]*r[i];

                       sqrt_normr=sqrt(normr);

         }

       while(sqrt_normr>eps);

             printf("U=\n");

       for(i=1;i<=p*p;i++)

           printf("%lf\n",U[i]);

             printf("sqrt_normr=%lf\n",sqrt_normr);

       printf("k=%d\n",k);

       system("pause");

       return 0;

}

SOR方法：

#include<stdio.h>

#include<stdlib.h>

#include<math.h>

int main()

{

      int  i,j,k;

      int  n,m,p;

      double  a,h,eps;

      double  w,sum,norm,sumB,sumr;

      double  u[400]={0},r[400]={0},B[400]={0};

      double  u1[400]={0},u0[400]={0},b[400]={0};

      double  A[400][400]={0};

      printf("n= \n");

      scanf("%d",&n);

      printf("eps= w= \n");

      scanf("%lf",&eps);

      scanf("%lf",&w);

      printf("eps=%lf\n",eps);

      printf("w=%lf\n",w);

      h=1.0/n;

      m=n*n;

      p=n-1;

      printf("%lf\n",h);

      a=h*h/4;

      printf("%d,%d,%lf,%lf\n",m,p,h,a);

          j=1;

      b[1]=a*sin(1.0/m)+2*a;

      b[n-1]=a*sin(1.0*p/m)+a*(p*p)+a*(1+n*n);

      for(i=2;i<=n-2;i++)

          b[i]=a*sin(1.0*i*j/m)+a*i*i;

      for(j=2;j<=n-2;j++)

      {

          b[(j-1)*p+1]=a*sin(1.0*j/m)+a*j*j;

          b[j*p]=a*sin(1.0*p*j/m)+a*(j*j+n*n);

          for(i=(j-1)*p+2;i<=j*p-1;i++)

              b[i]=a*sin((i-(j-1)*p)*j*1.0/m);

      }

      j=n-1;

      b[(p-1)*p+1]=a*sin(1.0*j/m)+p*p*a+a*(1+n*n);

    

      for(i=(p-1)*p+2;i<=p*p-1;i++)

          b[i]=a*sin(1.0*(i-(p-1)*p)*j/m)+1.0/4+a*(i-(p-1)*p)*(i-(p-1)*p);

      b[p*p]=a*sin(1.0*p*p/m)+(n*n+p*p)*2*a;

      A[1][1]=1+a;

      A[1][2]=-0.25;

      A[p*p][p*p-1]=-0.25;

      A[p*p][p*p]=1+a;

      for(i=2;i<=p*p-1;i++)

      {   A[i][i-1]=-0.25;

          A[i][i]=1+a;

          A[i][i+1]=-0.25;

      }

      for(i=n;i<=p*p;i++)

      {    A[i-n+1][i]=-0.25;

           A[i][i-n+1]=-0.25;

      }

      for(j=2;j<=n-1;j++)

      {   A[(j-1)*p][(j-1)*p+1]=0.0;

          A[(j-1)*p+1][(j-1)*p]=0.0;

      }

              k=0;

      do

      {

          k=k+1;

          for(i=1;i<=p*p;i++)           //求B=-U*u0+b U用A代替

           {  sumB=0.0;

              for(j=i+1;j<=p*p;j++)

                  sumB=sumB+A[i][j]*u0[j];

              B[i]=-sumB+b[i];

           }

                    u1[1]=1.0*B[1]/A[1][1];

           for(i=2;i<=p*p;i++)   //求u1

            {  sum=0.0;

               for(j=1;j<=i-1;j++)

                   sum=sum+A[i][j]*u1[j];

               u1[i]=1.0*(B[i]-sum)/A[i][i];

            }

            for(i=1;i<=p*p;i++)        //求解u

                u[i]=w*u1[i]+(1-w)*u0[i];

            for(i=1;i<=p*p;i++)        //求u0

                u0[i]=u1[i];

            for(i=1;i<=p*p;i++)      //求残量r

             {   sumr=0.0;

                 for(j=1;j<=p*p;j++)

                     sumr=sumr+A[i][j]*u[j];

                 r[i]=sumr-b[i];

             }

            norm=0.0;  

            for(i=1;i<=p*p;i++)//求r的范数

                norm=norm+r[i]*r[i];

            norm=sqrt(norm);

      }

      while( norm>eps);

      printf("k=%d\n",k);

      printf("norm=%lf\n",norm);

      for(i=1;i<=p*p;i++)

          printf("%lf\n",u[i]);

      system("pause");

      return 0;

}

运行结果：                                  

参考文献

【1】   徐树方，高立，张平文，数值线性代数。北京大学出版社，2010.

【2】       周品，何正风，MATLAB数值分析。机械工业出版社，2009.

更多相关推荐：
大作业报告格式

供配电技术课程大作业报告书题目指导教师姓名学号日期机电工和系20xx20xx学年第2学期报告书格式要求一报告前置部分一摘要内容包括研究目的方法结果结论300字400字四部分二格式要求1中文摘要摘要黑体三号居中摘...

大作业报告格式

面向对象技术课程大作业设计报告书题目超市管理系统指导教师宋涛姓名学号日期20xx1122管理科学与工程学院20xx20xx学年第1学期一需求分析随着小型超市规模的发展不断扩大商品数量急剧增加商品的各种信息量也成...

作业报告

关于仙林大学城饮料消费行为调查分析报告调查背景仙林大学城现在拥有的院校主要有南京大学南京师范大学南京财经大学南京邮电大学南京中医药大学南京信息职业技术学院南京理工大学紫金学院应天学院南京森林警察学院南京体育学院...

大作业报告格式

数据库原理大作业作业题目小组成员指导教师提交时间一需求分析功能分析功能结构图提取数据过程以及结果二概念结构设计根据需求分析得到实体属性以及实体间联系用文字描述清楚最后得到ER图三逻辑结构设计根据ER图转换相应的...

毕业大作业实践报告

中国石油大学华东现代远程教育毕业大作业实践报告题目工商企业安全管理实践报告学习中心奥鹏学习中心年级专业网络09秋工商企业管理学生姓名学号实践单位南通腾宇混凝土有限公司实践起止时间20xx年4月20xx年7月中国...

大作业报告封面

下大作业报告学院名称东方学院专业计算机科学与技术题目学期13142班级学号姓名报告成绩答辩成绩教师姓名李秉璋20xx年6月程序设计基础1作业性质大作业程是学生在学习完程序设计基础下后为提高学生编写程序的能力要求...

毕业大作业(实践报告)

中国石油大学华东现代远程教育毕业大作业实践报告题目工商企业管理专业实践报告学习中心青岛校区学习中心年级专业工商企业管理学生姓名学号实践单位实践起止时间中国石油大学华东远程与继续教育学院完成时间年月日0中国石油大...

Web程序设计用大作业报告模板

武汉工业学院Web高级程序设计大作业报告专业学号090502227姓名日期一需求分析1引言随着网络与信息技术的发展很多陌生人之间都有了或多或少的联系如何更好地管理这些信息是没跟人必须面临的问题特别是那些很久没有...

作业报告 (1)

汕头大学移动编程导论作业报告APPInventor设计报告马克思刷题器游戏设计开课班号68427学号20xxxxxxxx姓名xxx年级专业xxxx完成时间xxxxxxx一任务与要求根据本学期学习的appinve...

web大作业报告

Web程序设计期末考核报告院系数学与计算机学院专业软件工程班级1202班学号120xx10308姓名袁琦指导老师蒋丽华20xx年5月16日基于ASPNET的电子通信录系统的设计与实现一需求分析为了掌握使用ADO...

java大作业报告

西安思源学院Java大作业实验名称贪吃蛇游戏设计班级11计算机科学与技术学号1111020xx1041111020xx120组员指导教师王振铎摘要3一题目411贪吃蛇游戏412目的413适应实际实践编程的能力4...

C++大作业个人报告

C大作业个人工作总结1个人分工1编写地图类CFixedMap和AGV类CAGV2编写类CPathPlan中将地图转换为有向图的成员函数MaptoGraph3编写类CPathPlan中的Putpath函数的文件输...

大作业报告（32篇）

热门关注

C++大作业报告

毕业大作业(实践报告)

毕业大作业(实践报告)

大作业报告

大作业报告 2

中国石油大学(华东)现代远程教育 毕业大作业(实践报告)

综合作业报告格式

20xx中国大学生就业报告发布

猜你喜欢

医疗管理会议记录

医师定期考核管理办法

电气工程大学生职业规划书

中央电大毕业生登记表 样本20xx年

关于珍惜生命的格言（31篇）

祖国在我心中主题演讲稿（32篇）