数值分析上机作业一(Jab SOR CG)
数值分析上机报告
数值分析上机实践报告目录第一章标准题部分 (3)1.1 第一题 (3)1.1.1解题的理论依据、算法进行分析及应用条件 (3)1.1.2 MATLAB计算程序 (4)1..1.3 MATLAB运行结果 (7)1..1.4 结果分析与讨论 (7)1.2 第二题 (7)1.2.1解题的理论依据、算法进行分析及应用条件 (7)1.2.2 MATLAB计算程序 (8)1.2.2 MATLAB运行结果 (9)1.2.3 结果分析与讨论 (9)1.3 第三题 (9)1.3.1解题的理论依据、算法进行分析及应用条件 (9)1.3.2 MATLAB计算程序 (10)1.3.3 MATLAB运行结果 (11)1.3.4 结果分析与讨论 (11)1.4 第四题 (11)1.4.1解题的理论依据、算法进行分析及应用条件 (11)1.4.2 MATLAB计算程序 (12)1.4.3 MATLAB运行结果 (13)1.4.4 结果分析与讨论 (13)第二章自主题部分 (14)钢筋混凝土偏心受压柱配筋设计 (14)2.1 引言 (14)2.2 问题研究 (15)2.3 案例分析 (16)2.4 结论 (18)第一章 标准题部分1.1 第一题用三次样条插值的三弯矩法求f(-0.02)和f(2.56)。
1.1.1 解题的理论依据、算法进行分析及应用条件为使问题一般化，本题的三次样条求解使用了对于任意分化的三弯矩插值法。
若记各节点间距n x x x h j j j ,,2,111⋅⋅⋅=-=--，根据理论分析，三次样条插值函数为：))(6()()(6(6)(6)()(112112111131131-------------+--+-+-=j j j j j j j j j j j j jj j j h x x h M y h x x h M y h x x M h x x M x S其中),(1j j x x x -∈。
而式中的诸j M 可由如下的矩阵方程来确定：⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡=⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎣⎡N N d d d M M M10101111122212μλμλμ这里⎪⎪⎪⎪⎪⎪⎪⎩⎪⎪⎪⎪⎪⎪⎪⎨⎧-⋅⋅⋅=+=-=+=--'==---+='--=------+---+-)1,,2,1(1)](1[6),,(6)(6)(61111111111110000100N j h h h h h h y y h y h d x x x f h y y h y y h h d y h y y h d j j jj j j j j j N N N NN N j j j j j j j j j j j μλμ对于上述矩阵方程，由于其系数矩阵为一三对角阵，故我们用追赶法（Thomas 算法）来求解。
数值分析上机作业(MATLAB)
将系数矩阵 A 分解为：A=L+U+D
Ax=b
⇔ (D + L +U)x = b ⇔ Dx = −(L + U )x + b ⇔ x = −D −1(L + U )x + D −1b x(k +1) = −D −1 (L + U ) x(k ) + D −1b
输入 A，b 和初始向量 x
迭代矩阵 BJ , BG
否
ρ(B) < 1?
按雅各比方法进行迭代
否
|| x (k+1) − x(k) ||< ε ?
按高斯-塞德尔法进行迭代
否
|| x(k+1) − x (k ) ||< ε ?
输出迭代结果
图 1 雅各布和高斯-赛德尔算法程序流程图
1.2 问题求解
按图 1 所示的程序流程，用 MATLAB 编写程序代码，具体见附录 1。解上述三个问题 如下
16
-0.72723528355328
0.80813484897616
0.25249261987171
17
-0.72729617968010
0.80805513082418
0.25253982509100
18
-0.72726173942623
0.80809395746552
0.25251408253388
0.80756312717373
8
-0.72715363032573
0.80789064377799
9
-0.72718652854079
高等数值分析上机作业
运行结果：
2
拟合多项式为 y 0.9204x 0.4296x3 （2）用切比雪夫多项式拟合 arcsin x 程序如下；A 是正规方程的系数矩 阵，可以直接由公式得到,设拟合多项式为
y c0 c1x c2 (2x2 1) c3(4x3 3x)
下面将拟合曲线画出 x=-1:0.01:1; y=asin(x); plot(x,y,'r') hold on z=0.9204*x+0.4296*x.^3; plot(x,z,'g') hold on legend('y=arcsin(x)','y=0.9204x+0.4296x^3') B=zeros(4,1);
each
i
0,1,..., N
1.
程序代码：
建立M文件 function [x,y]=Euler1(fun,y0,h) x=0.1:h:0.5; y(1)=y0; for n=1:4
k1=feval(fun,x(n),y(n)); y(n+1)=y(n)+h*k1;
6
k2=feval(fun,x(n+1),y(n+1)); y(n+1)=y(n)+h*(k1+k2)/2; end x=x'; y=y'; 主程序 >> fun=inline('y-x.^2+1'); >> [x,y]=Euler1(fun,0.5,0.1)
2 S ( x)
2 S ( x)
a
f x 在子集 a,b中的最佳平方逼近函数。
取 j x x j , j 0,1,, n, 就有 span1, x, x2,, xn。对于任意的 Sx ，
数值分析上机实验报告
数值分析上机实验理学院11级统计01班41108030125鲁庆实验报告一一．实验名称误差与误差估计二．实验目的掌握数值运算的误差估计方法三．数学原理 1.绝对误差(*)e x设某一量的准确值为x ，近似值为x*，则x*与x 之差叫做近似值x*的绝对误差(简称误差)，记为*(*)*e e x x x ==- 2.绝对误差限适当小的正数，使|(*)||*|*e x x x ε=-≤则称*ε为近似值 x * 的绝对误差限。
(有时用*x x ε*=±表示近似值x *的精度或准确值的所在范围。
3.相对误差(*)r e x绝对误差与准确值之比*(*)*(*),0r r e x x xe e x x x x-===≠称为x *的相对 误差。
4.相对误差限(*)r x ε若指定一个适当小的正数 (*)r x ε，使|(*)||(*)|(*)||r r e x e x x x ε=≤则称(*)r x ε为近似值 x *的相对误差限。
5.有效数字若近似值x*的绝对误差限是某一位的半个单位，该位到x*的第一位非零数字一共有n 位，则称近似值x*有n 位有效数字，或说x*精确到该位。
6.绝对误差的运算：)()()(2121x x x x εεε+=± )()()(122121x x x x x x εεε+≈22122121+=x x x x x x x )()()(εεε (f(x))()(x)f x εε'≈四．实验内容1. 计算I n=e 1-⎰10nxe x 2dx (n=0,1,...)并估计误差。
解： >> I0 = exp(-1)*quad('(x.^0).*exp(x.^2)',0,1,10^(-10));>> vpa(I0,10) ans =.5380795069>> I1= exp(-1)*quad('(x.^1).*exp(x.^2)',0,1,10^(-10)); >> vpa(I1,10) ans =.3160602794>> I2 = exp(-1)*quad('(x.^2).*exp(x.^2)',0,1,10^(-10)); >> vpa(I2,10) ans =.2309602465>> I3 = exp(-1)*quad('(x.^3).*exp(x.^2)',0,1,10^(-10)); >> vpa(I3,10) ans =.1839397206>> I4 = exp(-1)*quad('(x.^4).*exp(x.^2)',0,1,10^(-10)); >> vpa(I4,10) ans =.1535596302>> I5 = exp(-1)*quad('(x.^5).*exp(x.^2)',0,1,10^(-10)); >> vpa(I5,10) ans =.1321205588>> I6 = exp(-1)*quad('(x.^6).*exp(x.^2)',0,1,10^(-10)); >> vpa(I6,10) ans =.1161009245>> I7 = exp(-1)*quad('(x.^7).*exp(x.^2)',0,1,10^(-10)); >> vpa(I7,10) ans =.1036383235>> I8 = exp(-1)*quad('(x.^8).*exp(x.^2)',0,1,10^(-10)); >> vpa(I8,10) ans =.9364676413e-1>> I9 = exp(-1)*quad('(x.^9).*exp(x.^2)',0,1,10^(-10)); >> vpa(I9,10) ans =.8544670595e-1 2.计算x255的值。
数值分析上机答案computer中科院
源程序如下：function[x_jac,x_gauss,x_sor,A,b]=diedai(n,w,e,N,Ain,bin) if nargin<3e=1e-5;N=1000;endif nargin<4N=1000;endif nargin<2error('²ÎÊý²»¹»');endA=hilb(n);b=A*ones(n,1);if nargin>4if nargin==5error('ÇëÊäÈëb');elseA=Ain;b=bin;endend%%Ñſ˱ȵü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ediedai_n=diedai_n+1;if diedai_n>Ndisp('Ñſ˱ȵü´úδ´ïµ½¾«¶È');break;endx(:,1)=x(:,2);for i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),1)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i);endendendx_jac.x=x(:,2);x_jac.n=diedai_n-1;x_jac.error=max(abs(b-A*x(:,2)));%%¸ß˹µü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('¸ß˹¡ªÈüµÂ¶ûµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,(i+1):n)*x((i+ 1):n,1))/A(i,i);elsex(i,2)=(b(i)-A(i,(i+1):n)*x((i+1):n,1))/A(i,i); endendendx_gauss.x=x(:,2);x_gauss.n=diedai_n-1;x_gauss.error=max(abs(b-A*x(:,2)));%%SORµü´údiedai_n=0;x=zeros(n,2);while max(abs(b-A*x(:,2)))>ex(:,1)=x(:,2);diedai_n=diedai_n+1;if diedai_n>Ndisp('SORµü´úδ´ïµ½¾«¶È');break;endfor i=1:nif i>1x(i,2)=x(i,1)+w*(b(i)-A(i,1:(i-1))*x(1:(i-1),2)-A(i,i:n)* x(i:n,1))/A(i,i);elsex(i,2)=x(i,1)+w*(b(i)-A(i,i:n)*x(i:n,1))/A(i,i);endendendx_sor.x=x(:,2);x_sor.n=diedai_n-1;x_sor.error=max(abs(b-A*x(:,2)));end运行结果：由于SOR的松弛因子取1时与高斯-赛德尔迭代相同，故在矩阵n取8,10时省去了SOR迭代松弛因子取1的步骤。
东南大学出版社第二版《数值分析》上机作业答案(前三章)
for (i=k+1;i<N;i++) // { lik=a[i][k]/a[k][k]; //实施消去过程，得到上三角系数增广矩阵 for (j=k;j<M;j++) // { a[i][j]=a[i][j]‐lik*a[k][j]; // } } } cout<<"经列主元高斯消去法得到的数组为:"<<endl; // for (b=0;b<N;b++) // { cout<<endl; //输出经过列主元消去法处理过的系数增广矩阵 for (c=0;c<M;c++) { cout<<setw(7)<<a[b][c]; // } } cout<<endl; double x[N]; // double s; int f,g; x[N‐1]=a[N‐1][M‐1]/a[N‐1][N‐1]; // for (f=N‐2;f>=0;f‐‐) // { s=0; for (g=f+1;g<N;g++) //由上三角形的系数增广矩阵求出方程组的解 { s=s+a[f][g]*x[g]; // } x[f]=(a[f][N]‐s)/a[f][f]; // } cout<<"方程组的解为:"<<endl; for (b=0;b<N;b++) //输出方程组的解 {
1
当 n=10000 时，s3=0.7499 Press any key to continue （分析 S1 的 6 位数字中，有效位数为 4 位； S2 的所有数字都是有效数字。 ） 当 n=1000000 时，s1=‐14.2546 当 n=1000000 时，s2=‐14.2551 当 n=1000000 时，s3=0.749999 Press any key to continue （分析： S1 的 6 位数字中，没有有效数字； S2 的 6 位数字中，没有有效数字。 ） 由运行结果可知，当精度比较低时，按从大数开始累加到小数的计算结果的精度低于按从小数 累加到大数的计算结果的精度。 至于当 n=1000000 时，S1 和 S2 得出了负数结果，可能是由于循环次数过多，导致数据溢出， 从而得出错误结果。 习题 2 20．程序如下： //给定误差限为：0.5e‐6 //经过试算得当 delta 最大取道 0.7745966 时，迭代得到的根都收敛于 0 #include <iostream.h> #include <math.h> void main () { double x,u; int count=0; u=10.0; cout<<"请输入 x 的初值"<<endl; cin>>x; for (count=0;abs(u)>5;count++) { x=x‐(x*x*x‐3*x)/(3*(x*x‐1)); u=10000000*x; if(count>5000) { cout<<"迭代结果不收敛于 0！"<<endl; break; } } cout<<"x="<<x<<endl<<endl;
数值分析上机实验
目录1 绪论 (1)2 实验题目(一) (2)2.1 题目要求 (2)2.2 NEWTON插值多项式 (3)2.3 数据分析 (4)2.3.1 NEWTON插值多项式数据分析 (4)2.3.2 NEWTON插值多项式数据分析 (6)2.4 问答题 (6)2.5 总结 (7)3 实验题目(二) (8)3.1 题目要求 (8)3.2 高斯-塞德尔迭代法 (8)3.3 高斯-塞德尔改进法—松弛法 (9)3.4 松弛法的程序设计与分析 (9)3.4.1 算法实现 (9)3.4.2 运算结果 (9)3.4.3 数据分析 (11)4 实验题目(三) (13)4.1 题目要求 (13)4.2 RUNGE-KUTTA 4阶算法 (13)4.3 RUNGE-KUTTA 4阶算法运算结果及数值分析 (14)总结 (16)附录A (17)1绪论数值分析是计算数学的一个主要部分，它主要研究各类数学问题的数值解法，以及分析所用数值解法在理论上的合理性。
实际工程中的数学问题非常复杂，所以往往需要借助计算机进行计算。
运用数值分析解决问题的过程：分析实际问题，构建数学模型，运用数值计算方法，进行程序设计，最后上机计算求出结果。
数值分析这门学科具有面向计算机、可靠的理论分析、好的计算复杂性、数值实验、对算法进行误差分析等特点。
本学期开设了数值分析课程，该课程讲授了数值分析绪论、非线性方程的求解、线性方程组的直接接法、线性方程组的迭代法、插值法、函数逼近与曲线拟合、数值积分和数值微分、常微分方程初值问题的数值解法等内容。
其为我们解决实际数学问题提供了理论基础，同时我们也发现课程中很多问题的求解必须借助计算机运算，人工计算量太大甚至无法操作。
所以学好数值分析的关键是要加强上机操作，即利用计算机程序语言实现数值分析的算法。
本报告就是基于此目的完成的。
本上机实验是通过用计算机来解答数值分析问题的过程，所用的计算工具是比较成熟的数学软件MATLAB。
大连理工大学矩阵与数值分析上机作业Word版
传播优3^ Word版文档•希望对您有帮助.可双击去除!矩阵与数值分析上机作业学校：大连理工大学学院：班级：姓名：学号：授课老师：注：编程语言Mat lab1.考虑计算给定向量的范数：输入向量x = (ri,x2,•■- ,-r n)r»输出胡I” ||创2，||広||oc・请编制一个通用程序，并用你编制的程序计算如下向量的范数：对“ =10, 100, 1000甚至更大的八计算其范数，你会发现什么结果？你能否修改你的程序使得计算结果相对精确呢？程序：Norm, m函数function s二Norm(x,m)给求向量x的范数師取1,2, in f分别表示仁2,无穷范数n=length(x);s=0;switch mcase 1 績-范数for i=1:ns二s+abs(x(i));endcase 2 %2-范数for i=1:ns 二s+x(i 厂2;ends=sqrt (s);case inf %无穷-范数s二max(abs(x));end计算向量X, y的范数Testi.mclear all;clc;n1=10；n2=100；n3=1000;x1=1./[1:n1],;x2=1./[1：n2]t ;x3=1./[1:n3]f;y1 = [1:n1]t;y2=[1:n2]-;y3=[1:n3],;disp(*n=10 时J ;dispC x 的 1 一范数：*) ;disp(Norm(x1, 1)); dispC x 的 2-范数：*) ;disp(Norm(x1,2)); dispC x 的无穷-范数:*) ;disp(Norm(x1, inf)); dispCy 的1-范数：’);disp(Norm(y1, 1)); dispC'y 的 2-范数:*) ;disp(Norm(y1,2));disp(* y 的无穷-范数：*) ;disp(Norm(y1, inf)); disp(*n=100 时J;dispC x 的1-范数：*) ;disp(Norm(x2, 1));dispCx 的 2-范数:*) ;disp(Norm(x2, 2)); disp(' x 的无穷-范數:*) ;disp(Norm(x2, inf)); disp(' y 的 1 一范数：*) ;disp(Norm(y2, 1)); dispC y 的 2-范数：*) ;disp(Norm(y2, 2));disp(* y 的无穷-范数:*) ;disp(Norm(y2, inf)); dispCn=1000 时’);disp(' x 的 1-范数：*) ;disp(Norm(x3, 1)); dispC x 的 2-范数：*) ;disp(Norm(x3, 2));dispC x 的无穷-范數:*) ;disp(Norm(x3, inf)); dispC y 的 1-范数：*) ;disp(Norm(y3, 1));dispC y 的 2-范数:*) ;disp(Norm(y3, 2)); dispC y 的无穷-范数:*) ;disp(Norm(y3, inf));y 的1-范数:500500： y 的2-范数：1・8272+004； y 的无穷-范数:10002. 考虑y = /(x)=芈也，其中定义/(O) = b 此时/仗)是连续函数.用此＞公式计算 当10-込10-】5］时的函数值，画出图像。
数值分析第一次上机
数值分析第一次上机姓名：学号：班级：一、已知数据表为：用四次拉格朗日插值多项式求f（0.596）的近似值。
二、已知数据表为：三、用于计算一、二大题的代码#include <iostream>#include <cmath>#include <fstream>using namespace std;class point{private:int n;double *x, *y, t, z; //x为存放结点的值；y为存放结点的函数值；t为插值点；z为所求结果public:point (int nn) //constructor{n = nn;x = new double[n]; //分配内存y = new double[n];}~point () //destructor{delete [] x, y;}void input (int count); //由文件读入count个数据点(x, y)double interp (double num,int count); //执行Pointange插值void output (); //输出插值点t处的近似值z到文件};void point::input (int count) //由文件读入count个数据点(x, y){int i;char str1[20];cout <<"输入文件名: ";cin >>str1;ifstream fin(str1);if (!fin){cout <<"无法打开该文件" <<str1 <<endl;exit (1);}for (i=0;i<count; i++) //读入数据(x,y){fin >>x[i];fin >>y[i];}fin.close ();}double point::interp (double num,int count) //执行Pointange插值{int i,j;double l;t=num;z=0;if(count <0){cout <<"error!"<<endl;return -1;}for (i=0;i<count;i++){l=1;for(j=0;j<count;j++){if(i!=j) l=l*(t-x[j])/(x[i]-x[j]);}z=z+l*y[i];}return z;}void point::output () //输出插值点t处的近似值z到文件{char str2[20];cout <<"输出文件名: ";cin >>str2;ofstream fout (str2,ios::app);if (!fout){cout <<"不能打开这个文件" <<str2 <<endl;exit(1);}fout <<endl <<t <<" " <<z <<endl;fout.close ();}void main () //main函数{point solution(10);double num;int count;cout <<"请输入插值结点的值!"<<endl;cin >>num;cout <<"请输入插值次数!"<<endl;cin >>count;solution.input (count+1); //由文件读入n个数据点(x, y)cout <<solution.interp(num,count+1)<<endl; //执行Pointange插值solution.output (); //输出插值点t处的近似值z到文}四、计算一、二大题的运行结果第一题：请输入插值结点的值!0.596请输入插值次数!4输入文件名: point11.txt0.631918输出文件名: point111.txtPress any key to continue其中，point11.txt的内容如下所示：0.4 0.410750.55 0.578150.65 0.696750.8 0.888110.9 1.02652point.111.txt的内容如下所示：0.596 0.631918第二题：请输入插值结点的值!11.5请输入插值次数!2输入文件名: point12.txt0.199369输出文件名: point111.txtPress any key to continue其中，point12.txt：11 0.19080912 0.207912130.224951point.112.txt：11.5 0.199369。
《数值分析》上机习题
《数值分析》上机习题1.用Newton法求方程X7 - 28X4 + 14=0在（0.1，1.9）中的近似根（初始近似值取为区间端点，迭代6次或误差小于0.00001）。
#include<stdio.h>#include<math.h>int main(){ float x1,x,f1,f2;static int count=0;x1=0.1 ;//定义初始值do{x=x1;f1=x*(x*x*x*x*x*x-28*x*x*x)+14;f2=7*x*x*x*x*x*x-4*28*x*x*x;//对函数f1求导x1=x-f1/f2; count++; }while(fabs(x1-x)>=1e-5);printf("%8.7f\n",x1); printf("%d\n",count);return 0;}2.已知函数值如下表：试用三次样条插值求f(4.563)及f’(4.563)的近似值。
include <stdio.h>#include <math.h>#define N 11main(){double B[N+1][N+1],m,x,u[N+1],y[N+1],c[N+1],d[N+1];double e[N+1]={2,0,4.15888308,6.5916738,8.3177664,9.6566268,10.750557,11.675460 6,12.47667,13.1833476,13.8155106,14.0155106};int i;x=4.563;B[0][0]=-2;B[0][1]=-4;B[N][N-1]=4;B[N][N]=2;for(i=1;i<N;i++){B[i][i-1]=1;B[i][i]=4;B[i][i+1]=1;}u[0]=B[0][0];y[0]=e[0];for(i=1;i<N;i++){m=B[i][i-1]/u[i-1];u[i]=B[i][i]-m*B[i-1][i];y[i]=e[i]-m*y[i-1];}c[N]=y[N]/u[N];for(i=N-1;i>=0;i--)c[i]=(y[i]-B[i][i+1]*c[i+1])/u[i];for(i=0;i<12;i++){m=fabs(x-i);if(m>=2)d[i]=0;else if(m<=1)d[i]=0.5*fabs(pow(m,3))-m*m+2.0/3;elsed[i]=(-1.0/6.0)*fabs(pow(m,3))+m*m-2*fabs(m)+4/3.0;} m=0;for(i=0;i<12;i++) m=m+c[i]*d[i];printf("f(%4.3f)=%f\n",x,m);printf("f'(4.563)=%lf\n",(c[4]-c[2])/2); }3.用Romberg 算法求)00001.0(sin )75(32314.1=+⎰ε允许误差dx x x x x . #include "stdafx.h" #include<stdio.h> #include<math.h> float f(float x) {float f=0.0;f=pow(3.0,x)*pow(x,1.4)*(5*x+7)*sin(x*x); return (f); } main() {int i=1,j,k,n=12;float T[12],a=1.0,b=3.0,s=0.0; T[0]=0.5*(b-a)*(f(a)+f(b));for(j=1;j<n-1;j++){ for(k=1;k<=pow(2,j-1);k++) s+=f(a+(2*k-1)*(b-a)/pow(2,j)); T[j]=0.5*(T[j-1]+(b-a)*s/pow(2,j-1)); s=0.0; }T[11]=(4*T[1]-T[0])/(float)3;for(;fabs(T[11]-T[0])>0.00001;i++) {T[0]=T[11];for(j=1;j<n-1-i;j++) T[j]=(pow(4,i)*T[j+1]-T[j])/(pow(4,i)-1); T[11]=(pow(4,i+1)*T[1]-T[0])/(pow(4,i+1)-1); }printf("%f\n",T[11]); }4. 用定步长四阶Runge-Kutta 求解一、程序要求用定步长四阶法求解 y1’=1y2’=y3 y3’=1000-1000y2-100y3(y1(0)=y2(0)=y3(0)=0) h=0.0005,打印yi(0.025),yi(0.045),yi(0.085),yi(0.1),(i=1,2,3)h =0.0005，打印y i (0.025) ，⎪⎪⎪⎪⎩⎪⎪⎪⎪⎨⎧===--===0)0(0)0(0)0(10010001000//1/321323321y y y y y dt dy ydt dy dt dyy i(0.045) ，y i(0.085) ，y i(0.1) ，(i=1，2，3)#include "stdafx.h"#include <stdio.h>#include <math.h>double F(double x,double y[4],double f[4]){f[1]=0*x+0*y[1]+0*y[2]+0*y[3]+1;f[2]=0*x+0*y[1]+0*y[2]+1*y[3]+0;f[3]=0*x+0*y[1]-1000*y[2]-100*y[3]+1000;return(1);}void main(){double F(double x,double y[4],double f[4]);doubleh=0.0005,x=0,Y[4],k[5][4],s[4],f[4],det,m[4]={0.025,0.045,0.085,0.1};int i,j,t; for(t=0;t<=3;t++)/*龙格-库塔算法*/{for(j=0;j<=3;j++)Y[j]=0; //每求一组值后将初值条件还原为0 for(i=1;i<=int(m[t]/h);i++){for(j=1;j<=3;j++)s[j]=Y[j];det=F(x,s,f);for(j=1;j<=3;j++)k[1][j]=h*f[j]; /*四阶古典公式中的k值和求和的计算*/ for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[1][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[2][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+0.5*k[2][j];det=F(x+0.5*h,s,f);for(j=1;j<=3;j++)k[3][j]=h*f[j];for(j=1;j<=3;j++)s[j]=Y[j]+k[3][j];det=F(x+h,s,f);for(j=1;j<=3;j++)k[4][j]=h*f[j];for(j=1;j<=3;j++)Y[j]=Y[j]+(k[1][j]+2*k[2][j]+2*k[3][j]+k[4][j])/6;x+=h;} for(j=1;j<=3;j++)printf("y[%d](%f)=%f ",j,m[t],Y[j]); printf("\n"); } } 5.⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎣⎡=40.00001 4.446782 2.213474- 0.238417 1.784317 0.037585- 1.010103- 3.123124 2.031743- 4.446782 30.719334 3.123789 1.103456- 2.121314 0.71828- 0.336993 1.112348 3.067813 2.213474- 3.123789 14.7138465 0.103458- 3.111223- 2.101023 1.056781- 0.784165- 1.7423820.238417 1.103456- 0.103458- 9.789365 0.431637 3.741856- 1.836742 1.563849 0.718719 1.7843172.1213143.111223- 0.431637 19.8979184.101011 2.031454 2.189736 0.113584-0.037585- 0.71828- 2.101023 3.741856- 4.101011 27.108437 3.123848 1.012345- 1.112336 1.010103- 0.336993 1.056781- 1.836742 2.031454 3.123848 15.567914 3.125432- 1.061074- 3.123124 1.112348 0.784165- 1.563849 2.189736 1.012345- 3.125432- 19.141823 2.115237 2.031743- 3.067813 1.742382 0.718719 0.113584- 1.112336 1.061074- 2.115237 12.38412A Tb )5.6784392- 4.719345 1.1101230 86.612343- 1.784317 0.84671695 25.173417- 33.992318 2.1874369(= 用列主元消去法求解Ax=b 。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
数值分析上机作业一 数值求解正方形域上的Poisson方程边值问题
解：由椭圆型第一边值问题的五点差分格式得到线性方程组为 2,1,1,,1,10,1,,0,141,(1)(1),(1)(1)01 , ijijijijijijjNjiiiiNijuuuuuhfijNuuyyjhjhuuxxihihijN
， 写成矩阵形式Au=f，其中
11,12,1,121,22,2,21,2,,(,,,),(,,,),,(,,,)TTTNNNNNNNvuuuvuuuvuuu 211,12,13,11,1,10,11,02,03,01,0,01,1221,22,2,20,21,221,2,3,1,,0,1,12,13,11,1(,,,,)(,,,,,)(,,,)(,0,,)(,,,,,)(,,,,,TTNNNNNTTNNTNNNNNNNNNNNNNNbhfffffuuuuuuubhfffuubhfffffuuuuuu
,11,)TNNNNu
,1,999,0.10.01 1,,1,2,...,1ijhNhfijNN取或则或 下面编写求解线性方程组Au=f的算法程序，用下列方法编程计算，并比较计算速度。 1) 用Jacobi迭代法求解线性方程组Au=f。 2) 用SOR迭代法求解线性方程组Au=f,用试算法确定最佳松弛因子。 3) 用共轭斜量法求解差分方程组Au=f。 1）Jacobi迭代法： function [U,k,er,t]=Jcb(n) % Jacobi迭代法 % U表示方程组的解；h表示步长；A表示迭代矩阵；k表示迭代次数；n表示非边界点数 % f表示线性方程组A*U=f的右端矩阵f；e表示允许误差界；er表示迭代误差 % t表示计算时间 h=1/(n+1); f(2:n+1,2:n+1)= h^2; %初始化f U=zeros(n+2); %初始化U矩阵，给出边界条件
for p=1:n+2 %由于matlab中下标规定从1开始，所以公式中所有对应下标都加1
U(1,p)=(p-1)*h*(1-(p-1)*h); U(n+2,p)=(p-1)*h*(1-(p-1)*h); U(p,1)=(p-1)*h*(1-(p-1)*h);
2222(,)1,0,1(0,)(1,)(1),01(,0)(,1)(1),01uufxyxyxyuyuyyyyuxuxxxx
1122
NN
vbvbufvb
4114114iiNNA
11
22
NN
AIIAAIIA
U(p,n+2)=(p-1)*h*(1-(p-1)*h); end A=U; e=0.000000001; %设置误差界
tic; for k=1:10000 %迭代求解
er=0; for i=2:n+1 for j=2:n+1 A(i,j)=(U(i-1,j)+U(i+1,j)+U(i,j+1)+U(i,j-1)+f(i,j))/4; er=er+abs(A(i,j)-U(i,j)); %估计当前误差 end end U=A; %将矩阵A的值赋予U er=er/n^2; if erend end t=toc; %获得计算时间
end
>> [U,k,er,t]=Jcb(9)
k = 322 er = 9.9052e-10 t = 0.0033 2）SOR function [U,k,er,t]=SOR(n,w) % SOR迭代法 % U 表示方程组的解；h 表示步长； k表示迭代次数；n 表示非边界点数 % f 表示线性方程组A*U=f的右端矩阵f ；e 表示允许误差界；er 表示迭代误差 % t 表示计算时间 h=1/(n+1); f(2:n+1,2:n+1)= h^2; %初始化f U=zeros(n+2); %初始化U矩阵，给出边界条件
for p=1:n+2 %由于matlab中下标规定从1开始，所以公式中所有对应下标都加1
U(1,p)=(p-1)*h*(1-(p-1)*h); U(n+2,p)=(p-1)*h*(1-(p-1)*h); U(p,1)=(p-1)*h*(1-(p-1)*h); U(p,n+2)=(p-1)*h*(1-(p-1)*h); end e=0.000000001; %设置误差界
tic; for k=1:10000 %迭代求解
er=0; for i=2:n+1 for j=2:n+1 Ub=U(i,j); U(i,j)=w*((U(i-1,j)+U(i+1,j)+U(i,j+1)+U(i,j-1)+f(i,j))/4)+(1-w)*U(i,j); er=er+abs(Ub-U(i,j)); %估计当前误差
end end er=er/n^2; if er%判断是否达到计算精度，如果达到则退出循环
end end t=toc; %获得计算时间
end
>> [U,k,er,t]=SOR(9,1.5)
k = 46 er = 8.8088e-10 t = 9.3534e-04 最佳松弛因子： function [goodw,goodk]=zuijiaw(n) %计算最佳松弛w %精度为0.001 [U,k,er,t]=SOR(n,0.01); %调用超松弛迭代函数，先取一个w的初值，得到k
goodk=k;goodw=0.01; for w=0.01:0.001:2 %在0
[U,k,er,t]=SOR(n,w); if k<=goodk goodk=k; goodw=w; end end end
>> [goodw,goodk]=zuijiaw(9) goodw = 1.5320 goodk = 33 3）共轭斜量法 function A=xishu(N) %存储离散化后非边界点构成的系数矩阵 A=zeros(N^2); for i=1:N^2 A(i,i)=4; if i+NA(i,i+N)=-1; A(i+N,i)=A(i,i+N); end if mod(i,N)~=0 A(i,i+1)=-1; A(i+1,i)=A(i,i+1); end end end
function [U,k,er,t]=cg(n) %共轭斜量法求解线性方程组 % 初始输入变量 a为系数矩阵 b为方程右端向量 x为输入的初始解 k为迭代次数 er表示迭代误差 h=1/(n+1); U=zeros(n+2); %初始化U矩阵
for p=1:n+2 U(1,p)=(p-1)*h*(1-(p-1)*h); U(n+2,p)=(p-1)*h*(1-(p-1)*h); U(p,1)=(p-1)*h*(1-(p-1)*h); U(p,n+2)=(p-1)*h*(1-(p-1)*h); end b0=h^2; %初始化b b(1,1)=b0+U(2,1)+U(1,2); for i=2:n-1 b(i,1)=b0+U(i+1,1); end b(n,1)=b0+U(n+1,1)+U(n+2,2);
for j=2:n-1 b((j-1)*n+1,1)=b0+U(1,j+1); b((j-1)*n+2:j*n-1,1)=b0; b(j*n,1)=b0+U(n+2,j+1); end
b(n*(n-1)+1,1)=b0+U(1,n+1)+U(2,n+2); for i=n*(n-1)+2:n^2-1 b(i,1)=b0+U(i-n*(n-1)+1,n+2); end b(n^2,1)=b0+U(n+1,n+2)+U(n+2,n+1);
a=xishu(n); e=0.000000001; %设置误差界
tic; for i=1:n^2 x(i,1)=0; end r=b-a*x; p=r; er=r'*r; for k=1:10000 w=a*p; t=er/(p'*w); x=x+t*p; r=r-t*w; q=r'*r; s=q/er; p=r+s*p; er=q; if erend end t=toc;
xx=reshape(x,n,n); U(2:n+1,2:n+1)=xx; end