Matlab大数值计算题目

传热学作业数值计算1 程序内容： 数值计算matlab程序内容：>> tw1=10; % 赋初值赋初值 tw2=20; c=1.5; p2=20; p1=c*p2; L2=40; L1=c*L2; deltaX=L2/p2; a=p2+1; b=p1+1; =ones(a,b)*5; m1=ones(a,b); m1(a,2:b-1)=zeros(1,b-2); m1(2:a,1)=zeros(a-1,1); m1(2:a,b)=zeros(a-1,1); m1(1,:)=ones(1,b)*2; k=0; max1=1.0; tn= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw1+tw2*sin(pi*(j-1)/(b-1)); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 t2=ones(a,b); %求解析温度场求解析温度场for i=a:-1:1 for j=1:1:b y=deltaX*(a-i); x=deltaX*(j-1); t2(i,j)=tw1+tw2*sin(pi*x/L1)*(sinh(pi*y/L1))/(sinh(pi*L2/L1)); end end t2 迭代次数k =706 数值解温度场 数值解每次迭代的最大误差max1 =9.8531e-07 解析温度场 t2 解析温度场取第11行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差数值计算matlab 程序内容：程序内容： >> tw1=10; tw2=20; c=1.5; p2=20; p1=c*p2; L2=20; deltaX=L2/p2; L1=c*L2; a=p2+1; b=p1+1; =ones(a,b)*5; m1=ones(a,b); m1(a,2:b-1)=zeros(1,b-2); m1(2:a,1)=zeros(a-1,1); m1(2:a,b)=zeros(a-1,1); m1(1,:)=ones(1,b)*2; k=0; max1=1.0; tn= ; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n= (i,j); switch m case 0 tn(i,j)=tw1; case 1 tn(i,j)=0.25*(tn(i,j+1)+tn(i,j-1)+tn(i+1,j)+tn(i-1,j)); case 2 tn(i,j)=tw2; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k max1 tx=ones(a,b); for i=1:1:a for j=1:1:b y=(a-i)*deltaX; x=(j-1)*deltaX; m=sym('m'); g=(((-1)^(m+1)+1)/m)*sin(m*pi*x/L1)*sinh(m*pi*y/L1)/sinh(m*pi*L2/L1); h=symsum(g,m,1,100); tx(i,j)=2*h*(tw2-tw1)/pi+tw1; end end tx 迭代次数k = 695 数值解温度场 数值解每次迭代的最大误差max1 =9.8243e-07 解析温度场 tx = 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点行的解析解的直线，散点为其数值解的点 曲线为第11行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像： 数值温度场 图像图像 ：解析温度场tx图像：数值解与解析解的误差数值解与解析解的误差程序内容： 数值计算matlab程序内容：>> t0=90; =10; L=10; c=0.25; p2=20; p1=p2/c; B=c*L; d=0.5*B; h=10; a=p2+1; b=p1+1; deltaX=B/p2; lambda=160; Bi=h*deltaX/lambda; =ones(a,b)*10; m1=ones(a,b)*3; m1(2:a-1,1)=zeros(a-2,1); m1(a,2:b-1)=ones(1,b-2); m1(1,2:b-1)=ones(1,b-2)*6; m1(2:a-1,b)=ones(a-2,1)*2; m1(1,b)=ones(1,1)*4; m1(a,b)=ones(1,1)*5; m1(1,1)=7; m1(a,1)=8; tn= ; max1=1.0; k=0; while ( max1>1e-6) k=k+1; max1=0; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=t0; case 1 tn(i,j)=(2*tn(i-1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 2 tn(i,j)=(2*tn(i,j-1)+tn(i-1,j)+tn(i+1,j)-4* )/(4+2*Bi)+ ; case 3 tn(i,j)=0.25*(tn(i,j-1)+tn(i,j+1)+tn(i-1,j)+tn(i+1,j)); case 4 tn(i,j)=(tn(i,j-1)+tn(i+1,j)-2* )/(2*Bi+2)+ ; case 5 tn(i,j)=(tn(i,j-1)+tn(i-1,j)-2* )/(2*Bi+2)+ ; case 6 tn(i,j)=(2*tn(i+1,j)+tn(i,j-1)+tn(i,j+1)-4* )/(4+2*Bi)+ ; case 7 tn(i,j)=t0; case 8 tn(i,j)=t0; end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k ta=ones(a,b); Bi1=h*d/lambda; sbi=sqrt(Bi1); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end x=deltaX*(j-1); ta(i,j)=(cosh(sbi*(L-x)/d)+sbi*sinh(sbi*(L-x)/d))*(t0- )/(cosh(sbi*L/d)+sbi*sinh(sbi*L/d))+ ; end end ta 迭代次数k =1461 数值解温度场 解析温度场 ta 解析温度场行的解析解和数值解的点取第11行的解析解和数值解的点曲线为第11行的解析解的直线，散点为其数值解的点行的解析解的直线，散点为其数值解的点解析解 第11行的误差=[数值解(11行) –解析解(11行)]/解析解图像如下图像如下数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差程序内容：数值计算matlab程序内容：>> tw=10; L2=15; c=0.75; L1=L2/c; p2=24 ; p1=p2/c; deltaX=2*L2/p2; a=p2+1; b=p1+1; lambda=16; qv0=24; =ones(a,b)*5; m1=ones(a,b); m1(1,:)=zeros(1,b); m1(2:a,b)=zeros(a-1,1); m1(2:a,1)=zeros(a-1,1); m1(a,2:b-1)=zeros(1,b-2); tn= ; max1=1.0; k=0; while(max1>1e-6) max1=0; k=k+1; for i=1:1:a for j=1:1:b m=m1(i,j); n=tn(i,j); switch m case 0 tn(i,j)=tw; case 1 tn(i,j)=0.25*(tn(i-1,j)+tn(i+1,j)+tn(i,j-1)+tn(i,j+1)+qv0*(deltaX^2)/lambda); end er=abs(tn(i,j)-n); if er>max1 max1=er; end end end =tn; end k; tx=ones(a,b); for i=1:1:a for j=1:1:b if i>(a+1)/2 y=-(i-(a+1)/2)*deltaX; else y=((a+1)/2-i)*deltaX; end if j>(b+1)/2 x=(j-(b+1)/2)*deltaX; else x=-((b+1)/2-j)*deltaX; end m=sym('m'); xi=(2*m-1)*pi/2; g=((-1)^m)/(xi^3)*(cosh(xi*y/L1)/cosh(xi*L2/L1))*cos(xi*x/L1); h=symsum(g,m,1,100); tx(i,j)=2*qv0*L1^2/lambda*h+qv0*(L1^2-x^2)/(2*lambda)+tw; end end tx 数值温度场 解析温度场tx 取第13行的解析解和数值解的点行的解析解和数值解的点行的解析解的直线，散点为其数值解的点曲线为第13行的解析解的直线，散点为其数值解的点解析解 第13行的误差=[数值解(13行) –解析解(13行)]/解析解数值温度场图像数值温度场图像 解析温度场图像解析温度场图像数值解与解析解的误差数值解与解析解的误差。

佛山科学技术学院实验报告课程名称_______________ 数值分析________________________实验项目_______________ 数值积分____________________专业班级机械工程姓名余红杰学号2111505010 指导教师陈剑成绩日期月日一、实验目的b1、理解如何在计算机上使用数值方法计算定积分 a f ""X的近似值；2、学会复合梯形、复合Simpson和龙贝格求积分公式的编程与应用。

3、探索二重积分.11 f （x, y）dxdy在矩形区域D = ｛（ x, y） | a _ x _ b, c _ y _ d｝的数值D积分方法。

二、实验要求（1）按照题目要求完成实验内容；（2）写出相应的Matlab程序；（3）给出实验结果（可以用表格展示实验结果）；（4）分析和讨论实验结果并提出可能的优化实验。

（5）写出实验报告。

三、实验步骤1、用不同数值方法计算积xln xdx =-- 0 9（1）取不同的步长h，分别用复合梯形及复合辛普森求积计算积分，给出误差中关于h的函数，并与积分精确值比较两公式的精度。

（2）用龙贝格求积计算完成问题（1 ）。

2、给出一种求矩形区域上二重积分的复化求积方法，然后计算二重积分..e"y dxdy，其中积分区域D二｛0乞x岂1,0岂y乞1｝。

1.%lnt_t.m复化梯形：function F = Int_t(x1,x2,n)%复化梯形求积公式% x1,x2为积分起点和中点%分为n个区间，没选用步长可以防止区间数为非整数。

%样点矩阵及其函数值：x = lin space(x1,x2 ,n+1);y = f(x);m = len gth(x);%本题中用Matlab计算端点位置函数值为NaN,故化为零: y(1) = 0;y(m) = 0;%算岀区间长度，步长h：h = (x2 -x1)/n;a = [1 2*o nes(1,m-2) 1];%计算估计的积分值：F = h/2*sum(a.*y);%f.mfun cti on y = f(x)y = sqrt(x).*log(x);%run 11.mclc,clear;%分为10个区间，步长0.1的积分值：F = In t_t(0,1,10);F10 = F%分为100个区间F = In t_t(0,1,100);F100 = F%误差计算W10 = abs((-4/9)-F10);W100 = abs((-4/9)-F100);W = [W10 W100]%复化辛普森：%l nt_s.mfun cti on F = In t_s(x1,x2 ,n)%复化梯形求积公式% x1,x2区间，分为n个区间。

1、设⎥⎦⎤⎢⎣⎡++=)1(sin35.0cos2xxxy，把x=0~2π间分为101点，画出以x为横坐标，y为纵坐标的曲线。

第一题的matlab源程序：①考虑cos（x）为一个整体，然后乘以中括号里面的全部x=0:2*pi/100:2*pi; %x的步长以及范围从0到2*pi y=cos(x).*(0.5+3*sin(x)./(1+x.^2)); %y的表达式plot(x,y)%画出图形图如下：②考虑对整体求解cos，先求x乘以括号中的部分x=0:2*pi/100:2*pi; %x的步长以及范围从0到2*pi y=cos(x.*(0.5+3*sin(x)./(1+x.^2))); %y的表达式plot(x,y) %画出图形图如下：2、产生8×6阶的正态分布随机数矩阵R1, 求其各列的平均值和均方差。

并求该矩阵全体数的平均值和均方差。

第二题的matlab源程序如下：R1=randn(8,6) %产生正态分布随机矩阵R1 =1.0933 -0.7697 1.5442 -0.1924 1.4193 0.21571.1093 0.3714 0.0859 0.8886 0.2916 -1.1658-0.8637 -0.2256 -1.4916 -0.7648 0.1978 -1.14800.0774 1.1174 -0.7423 -1.4023 1.5877 0.1049-1.2141 -1.0891 -1.0616 -1.4224 -0.8045 0.7223-1.1135 0.0326 2.3505 0.4882 0.6966 2.5855-0.0068 0.5525 -0.6156 -0.1774 0.8351 -0.66691.5326 1.1006 0.7481 -0.1961 -0.2437 0.1873aver=(sum(R1(1:end,1:end)))./8 %产生各行的平均值aver =0.0768 0.1363 0.1022 -0.3473 0.4975 0.1044a=std(R1(1:end,1:end)) %产生各行的均方差也就是标准差a =1.0819 0.8093 1.3456 0.8233 0.8079 1.2150aver1=(sum(R1(:)))./48 %全体数的平均值aver1 =0.0950b=std(R1(:)) %全体数的均方差即标准差b =1.01033、设x=rcost+3t,y=rsint+3,分别令r=2,3,4，画出参数t=0~10区间生成的x~y 曲线。

1、设⎥⎦⎤⎢⎣⎡++=)1(sin35.0cos2xxxy，把x=0~2π间分为101点，画出以x为横坐标，y为纵坐标的曲线。

第一题的matlab源程序：①考虑cos（x）为一个整体，然后乘以中括号里面的全部x=0:2*pi/100:2*pi; %x的步长以及范围从0到2*pi y=cos(x).*(0.5+3*sin(x)./(1+x.^2)); %y的表达式plot(x,y)%画出图形图如下：②考虑对整体求解cos，先求x乘以括号中的部分x=0:2*pi/100:2*pi; %x的步长以及范围从0到2*pi y=cos(x.*(0.5+3*sin(x)./(1+x.^2))); %y的表达式plot(x,y) %画出图形图如下：2、产生8×6阶的正态分布随机数矩阵R1, 求其各列的平均值和均方差。

并求该矩阵全体数的平均值和均方差。

第二题的matlab源程序如下：R1=randn(8,6) %产生正态分布随机矩阵R1 =1.0933 -0.7697 1.5442 -0.1924 1.4193 0.21571.1093 0.3714 0.0859 0.8886 0.2916 -1.1658-0.8637 -0.2256 -1.4916 -0.7648 0.1978 -1.14800.0774 1.1174 -0.7423 -1.4023 1.5877 0.1049-1.2141 -1.0891 -1.0616 -1.4224 -0.8045 0.7223-1.1135 0.0326 2.3505 0.4882 0.6966 2.5855-0.0068 0.5525 -0.6156 -0.1774 0.8351 -0.66691.5326 1.1006 0.7481 -0.1961 -0.2437 0.1873aver=(sum(R1(1:end,1:end)))./8 %产生各行的平均值aver =0.0768 0.1363 0.1022 -0.3473 0.4975 0.1044a=std(R1(1:end,1:end)) %产生各行的均方差也就是标准差a =1.0819 0.8093 1.3456 0.8233 0.8079 1.2150aver1=(sum(R1(:)))./48 %全体数的平均值aver1 =0.0950b=std(R1(:)) %全体数的均方差即标准差b =1.01033、设x=rcost+3t,y=rsint+3,分别令r=2,3,4，画出参数t=0~10区间生成的x~y 曲线。

Matlab 在化工数值计算中的应用(提纲) 基础知识Command Window 指令窗简介最简单的计算器使用方法加减乘除和幂运算符、矩阵的输入形式、常见表达式形式数值、变量和表达式数值的表示方法(十进制、科学记数)、变量命名规则(对大小写敏感，变量名的第一个字母必须为英文字母，不得含空格，但可含下划线链接符)、Matlab 默认的预定义变量(ans/inf/i 或j/pi/NaN 等)、复数和复数矩阵(把复数作为一个整体处理、real(z),imag(z),abs(z)/模,angle(z)/相角)。

例1：已知/612334,12,2i z i z i z e π=+=+=，并计算123/z z z z =计算结果的图形表示。

例2：画出衰减振荡曲线/3sin3t y e t -=及它的包络线并计算/30t y e -=。

t 的取值范围是[]04π， t=0:pi/80:4*pi; %定义自变量取值数组y0=exp(-t/3); %计算与自变量相应的y0数组y=exp(-t/3).*sin(3*t);%计算与自变量相应的y 数组plot(t,y,'-r',t,y0,':b',t,-y0,':b') %用不同颜色，不同线条绘制曲线数值计算结果的显示格式format/format short, format long, format short e, format long e, 标点符号的使用指令窗的常用控制指令clc, clear,edit, help, exit/quit, typeM 角本文件的编写与运算路径的制定帮助系统数值数组及其运算数组及其运算是Matlab 的核心内容2.1 一维数组的创建与赋值(逐个元素输入法、冒号生成法)；2.2 二维数组的创建与复制(直接输入法)2.3 执行数组运算的常用函数三角函数、反三角函数、幂指对函数、复数函数(abs,angle,conj （共厄复数）,imag,real)2.4 数组运算与矩阵运算A.’, A’,S./B,s*inv(B),A.^n,A^n,A.*B,A*B,A./B,A/B,f(A)注意运算符的小黑点。

1.圆盘上有如图所示的二十个数，请找出哪4个相邻数之和最大，并指出它们的位置和数值。

（10分）2011841361015217319716811149125解答：%1.圆盘上有如图所示的二十个数，请找出哪4个相邻数之和最大，并指出它们的位置和数值。

（10分）A=[1 18 4 13 6 10 15 2 17 3 19 7 16 8 11 14 9 12 5 20];% 程序位置规定：从1开始顺时针方向计数；NumA=size(A); Num=NumA(1,2); sum(1)=A(1); for i=1:(Num-3)sum(i)=A(i)+A(i+1)+A(i+2)+A(i+3); endmaxresult=max(sum(:));%找出4个相邻数之和最大值 maxresult %4个相邻数之和最大值 Position=i %四个数起始位置 FourNumber=A(1,i:(i+3)) %四个数的值及顺序 运行结果：maxresult =50Position =17FourNumber =9 12 5 202.甲、乙、丙三人上街买糖果。

三人都买好后，甲对乙、丙说，我可以按你们现有的糖果数再送你们每人一份。

甲送给乙、丙后，乙也按甲、丙现有的糖果数，送给甲、丙每人各一份糖果。

丙也如此送了甲、乙各一份。

互相赠送后，每人恰好各有64颗糖果。

问甲、乙、丙原来各买了多少糖果？（10分）解答：%由代数关系构造矩阵 a=[1 -1 -1;0 2 0;0 0 2]; b=a([2 1 3],:); b=b(:,[2 1 3]); c=a([2 3 1],:); c=c(:,[2 3 1]); d=64*ones(3,1); result=a\(b\(c\d))运行结果：result=104 56 323.求n S a aa aaa aaa a =++++ 的值。

a 的值为1~9之间的一个整数，n S 中每一项aaa a 共有n 位。

课程设计课程名称：设计题目：学号：姓名：完成时间：题目一：非线性方程求根 一 摘要非线性方程的解析解通常很难给出，因此非线性方程的数值解就尤为重要。

本实验通过使用常用的求解方法二分法和Newton 法及改进的Newton 法处理几个题目，分析并总结不同方法处理问题的优缺点。

观察迭代次数，收敛速度及初值选取对迭代的影响。

用Newton 法计算下列方程（1） 310x x --= ， 初值分别为01x =，00.45x =，00.65x =； （2） 32943892940x x x +-+= 其三个根分别为1,3,98-。

当选择初值02x =时给出结果并分析现象,当6510ε-=⨯，迭代停止。

解：1)采用MATLAB 进行计算；首先定义了Newton 法：function kk=newton(f,df,x0,tol,N)% Newton Method （牛顿法）% The first parameter f is a external function with respect to viable x.（第一个参数也就是本题所用的函数f ）% The second parameter df is the first order diffential function of fx.（第二个参数也就是本体所用函数f 的导数方程df ） % x0 is initial iteration point(初值). % tol is the tolerance of the loop （精度）.% N is the maximum number of iterations （循环上限）. x=x0;f0=eval(f);df0=eval(df); n=0;disp(' [ n xn xn+1 fn+1 ]'); while n<=N x1=x0-f0/df0; x=x1; f1=eval(f); X=[n,x0,x1,f1]; disp(X);if abs(x0-x1)<tolfprintf('The procedure was successful.') kk=X; return else n=n+1; x0=x1;f0=f1;endendif n==N+1fprintf('the method failed after N iterations. '),kk=0;End我们把Newton法存为.m格式的文件；之后我们运行程序：clear;clc;syms xf=x^3-x-1;df=diff(f,x);x=newton(f,df,1,0.0001,50);x会得到一下结果[ n xn xn+1 fn+1 ]0 1.0000 1.5000 0.87501.0000 1.5000 1.0625 -0.86302.0000 1.0625 1.4940 0.8408到第50次迭代时候会出现该问题：47.0000 1.4898 1.0814 -0.816748.0000 1.0814 1.4898 0.816749.0000 1.4898 1.0814 -0.816750.0000 1.0814 1.4898 0.8167the method failed after N iterations.x =0；同样测试x0=0.45、0.65得不出结果，判断出初值离真值太远，所以我们采用牛顿下山法进行计算迭代：我们定义了其中的f函数和df函数，并且分别存为.m格式的文件，其代码如下：f:function y=f(x)y=x^3-x-1;df:function y=df(x)y=3*x^2-1;之后我们定义newton下山法同时也存为.m的程序：function [x,i]=downnewton(f,df,x0,tol)k=0;i=1;disp(' [ n xn xn+1 fn+1 ]'); while(k==0)fx=feval('f',x0);dfx=feval('df',x0);t=0;u=1;while(t==0)dx=-fx/dfx;x1=x0+u*dx;fx1=feval('f',x1);fx0=feval('f',x0);if(abs(fx1)>abs(fx0));u=u/2;elset=1;endendX=[i,x0,x1,fx1];disp(X);if(abs(fx1)<tol)k=1;elsex0=x1;i=i+1;endendx=x1;i=i;end之后带入x0=0.45;downnewton('f','df',0.45,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.4500 -0.4155 -0.65622.0000 -0.4155 -0.5857 -0.61523.0000 -0.5857 -0.5754 -0.61514.0000 -0.5754 -0.5782 -0.61515.0000 -0.5782 -0.5773 -0.61516.0000 -0.5773 -0.5774 -0.61517.0000 -0.5774 -0.5773 -0.61518.0000 -0.5773 -0.5774 -0.61519.0000 -0.5774 -0.5774 -0.615110.0000 -0.5774 -0.5774 -0.615111.0000 -0.5774 1.3131 -0.049012.0000 1.3131 1.3248 0.000513.0000 1.3248 1.3247 0.0000ans =1.3247带入x0=0.6;downnewton('f','df',0.6,10^(-6))[ n xn xn+1 fn+1 ]1.0000 0.6000 1.1406 -0.65662.0000 1.1406 1.3668 0.18663.0000 1.3668 1.3263 0.00674.0000 1.3263 1.3247 0.00005.0000 1.3247 1.3247 0.0000ans =1.3247带入x0=1;downnewton('f','df',1,10^(-6))[ n xn xn+1 fn+1 ]1.0000 1.0000 1.5000 0.87502.0000 1.5000 1.3478 0.10073.0000 1.3478 1.3252 0.00214.0000 1.3252 1.3247 0.0000ans =1.32472)同样采用Newton下山法：重新定义f、df:f:function y=f(x)y=x^3+94*x^2-389*x+294;df:function y=df(x)y=3*x^2+188*x-389;再带入初值x0=2;downnewton('f','df',2,5*10^(-6))[ n xn xn+1 fn+1 ]1 2 -98 0ans =-98得出x=-98；分析：先画出该函数的图像；x=(-100:.1:100);ezplot('x^3+94*x^2-389*x+294',[-100 100]) 得出该图像如图：-100-80-60-40-20020406080024681012141618x 105xx 3+94 x 2-389 x+294根据牛顿法的几何解释，在x0=2的点做切线，与y 相交，交点的横坐标值为x=-98则结束了该现象。