矩阵算法程序

4x1x2x3x41x4xxx11234 xx4xx13412x1x2x34x411.用直接法求解；

2.用Jacobi迭代法求解；

3.分别取0.75,1.0,1.25,1.5，用SOR方法求解．比较迭代结果（与精确解比较）。 解：

1、用直接法求解 算法：Gauss列主元消去法是在Gauss消去法中增加选主元的过程，即在第k步（k=1,2,3,…）

消元时，首先在第k列主对角元以下（含对角元）元素中挑选绝对值最大的数（即为列主元），并通过初等行变换，使得该数位于主对角线上，然后再继续消元。 程序： /*gauss.m*/

function x=gauss(A,b,n) a=[A,b];/*a为增广阵*/ /*消去过程*/ for k=1:n-1

/*选主元*/ c=0; for q=k:n

if abs(a(q,k))>c c=a(q,k); l=q; end end

/*如果主元为0，则矩阵A不可逆*/ if abs(c)<1e-10 disp('error'); pause; exit; end

/*如果1不等于k，则交换第1行和第k行*/ if l~=k

for q=k:n+1

temp=a(k,q); a(k,q)=a(l,q); a(l,q)=temp; end end

/*计算第k步的元素值*/ for i=k+1:n

for j=k+1:n+1

a(i,j)=a(i,j)-a(i,k)/a(k,k)*a(k,j); end end end

/*回带过程*/

x(n)=a(n,n+1)/a(n,n); for i=n-1:-1:1 s=0;

for j=i+1:n

s=s+a(i,j)*x(j); end

x(i)=(a(i,n+1)-s)/a(i,i); end

/*返回方程组的解*/

fprintf('方程组的解为：')

>> A=[-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4]’ >> b=[1 1 1 1]’ >>n=4

>>gauss(A,b,n) 计算结果： 方程组的解为 x =

-1 -1 -1 -1

2、Jacobi迭代法

算法：对于线性方程组Ax=b，如果A为非奇异方程，则可将A 分解为：A=D-L-U其中D

为对角阵，其元素为A的对角元素，L与U为A 的下三角阵和上三角阵。于是Ax=b

化为：xk1D1(LU)xkD1b，其中BJD1(LU)，fDb。

matlab程序：

function x=jacobi(A,b,x0) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1); B=D\\(L+U); f=D\\b; x=B*x0+f; n=1;

while norm(x-x0,2)>=1.0e-10 x0=x;

1

x=B*x0+f; n=n+1; end

fprintf('迭代次数为：') n

fprintf('\\n')

fprintf('方程组的解为：')

>> A=[-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4]’ >> b=[1 1 1 1]’ >> x0=[0 0 0 0]'

计算结果： 迭代次数为： n = 79

方程组的解为： ans =

-0.99999999986515 -0.99999999986515 -0.99999999986515 -0.99999999986515

3、SOR法

算法：对于线性方程组Ax=b，如果A为非奇异方程，则可将A 分解为：A=D-L-U其中D为对

角阵，其元素为A的对角元素，L与U为A 的下三角阵和上三角阵。SOR法是在Gauss-Seidel迭代法的基础上引入松弛因子w，于是Ax=b化为：xk1Lwxkf， 其中Lw(DwL)1[(1w)DwU]，f(DwL)b Matlab程序：

function x=sor(A,b,x0,w) D=diag(diag(A)); U=-triu(A,1); L=-tril(A,-1);

Q=(D-w*L)\\((1-w)*D+w*U); f=(D-w*L)\\b*w; x=Q*x0+f; n=1;

while norm(x-x0)>=1.0e-10 x0=x; x=Q*x0+f; n=n+1; end

fprintf('迭代次数为：')

1

n

fprintf('\\n')

fprintf('方程组的解为：')

>> A=[-4 1 1 1;1 -4 1 1;1 1 -4 1;1 1 1 -4]’ >> b=[1 1 1 1]’ >> x0=[0 0 0 0]'

>> sor(A,b,x0,0.75) >> sor(A,b,x0,1.0) >> sor(A,b,x0,1.25) >> sor(A,b,x0,1.5)

w=0.75时计算结果为： 迭代次数为： n = 74

方程组的解为： x =

-0.99999999988656 -0.999999999506 -0.99999999990292 -0.99999999991020

w=1.0时计算结果为： 迭代次数为 n = 42

方程组的解为 x =

-0.99999999992282 -0.99999999993294 -0.99999999994173 -0.99999999994937 w=1.25时计算结果为： 迭代次数为 n = 21

方程组的解为 x =

-1.00000000000666 -0.99999999999062 -1.00000000000561 -0.99999999999973 w=1.5时计算结果为： 迭代次数为

n = 39

方程组的解为 x =

-1.00000000000227 -1.00000000001380 -0.99999999998230 -1.00000000001033

方程组的精确解为x=[1 1 1 1]’,0二、三次样条插值

算法： 三次样条插值是一种特殊类型分段三次插值，它具有一致收敛，充分光

滑和插值条件简单的特点。本题条件为自然条件s'(0)0.8,s'(10)0.2通过迭代计算各三次多项式的系数进行插值。

程序：

/*以下程序计算各三次多项式的系数*/

function S=scfit(x,y,dx0,dxn) N=length(x)-1; H=diff(x); D=diff(y)./H; A=H(2:N-1);

B=2*(H(1:N-1)+H(2:N)); C=H(2:N); U=6*diff(D); B(1)=B(1)-H(1)/2; U(1)=U(1)-3*(D(1)-dx0); B(N-1)=B(N-1)-H(N)/2; U(N-1)=U(N-1)-3*(dxn-D(N)); for k=2:N-1

temp=A(k-1)/B(k-1); B(k)=B(k)-temp*C(k-1); U(k)=U(k)-temp*U(k-1); end

M(N)=U(N-1)/B(N-1); for k=N-2:-1:1

M(k+1)=(U(k)-C(k)*M(k+2))/B(k); end

M(1)=3*(D(1)-dx0)/H(1)-M(2)/2; M(N+1)=3*(dxn-D(N))/H(N)-M(N)/2; for k=0:N-1

S(k+1,1)=(M(k+2)-M(k+1))/(6*H(k+1));

S(k+1,2)=M(k+1)/2;

S(k+1,3)=D(k+1)-H(k+1)*(2*M(k+1)+M(k+2))/6; S(k+1,4)=D(k); end

计算结果： ans =

-0.00851403685004 -0.00148596314996 0.80000000000000 0 -0.0044578449 -0.02702807370007 0.77148596314996 0.79000000000000 -0.00365440535039 -0.04040174204975 0.70405614740014 1.53000000000000 -0.040924414854 -0.051395810093 0.612244724946 2.19000000000000 0.10735236194454 -0.17413842554654 0.38678606360200 2.71000000000000 -0.26848495862962 0.14791866028708 0.36056629834254 3.03000000000000 0.42658747257395 -0.65753621560179 -0.14905125697216 3.27000000000000 -0.2673166618 0.62222620212007 -0.18436127045388 2.000000000000 0.05487225409078 -0.18136859287848 0.259633878770 3.06000000000000 0.05837591530307 -0.01675183060615 0.05837591530307 3.19000000000000 另外也可以使用matlab中的csape函数进行插值，命令如下：

>> x=[0:1:10];

>> y=[0 0.79 1.53 2.19 2.71 3.03 3.27 2. 3.06 3.19 3.29]; >>pp=csape(x,y,'complet',[0.8 0.2])；

>> breaks,coefs,nploys,ncofs,dim]=unmkpp(pp)/*显示节点、系数等*/ 得到的三次样条插值曲线如下图所示：

三、数值积分

1/e算法：采用变量替换，令tex，则原式=01dt，由于积分下限为0，算法中可令lnta1/2^n0，n10(或更大)，b1，采用复合辛普生公式，h(ba)/M，e节点数为2M+1,M可根据精度要求调整。 程序：

function s=simpson()

a=1/(1024*1024); /*n=20*/ b=1/exp(1);

h=(b-a)/200000; /*M=100000*/ s1=0; s2=0;

for k=1:100000 x=a+h*(2*k-1); s1=s1+1/log(x); end

for k=1:99999 x=a+2*h*k; s2=s2+1/log(x); end

s=h*(log(a)+log(b)+4*s1+2*s2)/3; s=-s;

计算结果为： ans =

0.219392325352

当取n=30,M=1000000时，计算结果为： ans =

0.219385205918

所以可以认为取n=20，M=100000，已经达到精度要求，因此



exdx0.2193923253 52x1四、求根

function A=zero()

Y=inline('sin(x)-(x.^2)/2','x'); x=-4:0.01:4; y=sin(x)-(x.^2)/2;

plot(x,y,'r');hold on;plot (x,zeros(size(x)),'b'); xlabel('x');ylabel('y');hold off title('sin(x)-0.5*x^2');

for i=1:2

[t1,y1]=ginput(1);

[t(i),yy(i),exitflag]=fzero(Y,t1) end A(:,1)=t'; A(:,2)=yy';

计算结果： ans =

0.00000000000000 0.00000000000000 1.40441482409243 0 即：方程的根为：

x1=1.40441482409243 x2=0

五、解微分方程

u2t3u,u(0)1,1、线性二步法

function A=erbufa() h=0.05; u(1)=0;

u(2)=2/9*exp(-0.15)+2/3*(0.05-1/3); for n=1:19 t(n)=h*(n-1); t(n+1)=h*n;

0t1

t(n+2)=h*(n+1);

m=2/3*h*(2*t(n+1)-3*u(n+1));

(n+2)=1/(1+0.25*h)*(u(n+1)+0.25/6*t(n+2)+m-5/12*h*(2*t(n)-3*u(n))); end A=[t' u']; plot(t',u','r'); legend('线性二步法')

计算结果： r =

0 0.05000000000000 0.10000000000000 0.15000000000000 0.20000000000000 0.25000000000000 0.30000000000000 0.35000000000000 0.40000000000000 0.45000000000000 0.50000000000000 0.55000000000000 0.60000000000000 0.65000000000000 0.70000000000000 0.75000000000000 0.80000000000000 0.85000000000000 0.90000000000000 0.95000000000000 1.00000000000000

0

0.00237955031668 0.00952256324446 0.01931097496112 0.031744316983 0.046693662634 0.0041178410 0.08367612820766 0.10549243349842 0.129381999495 0.15526910790694 0.18304099221326 0.21261683513033 0.24391296819485 0.27684960352693 0.31135065384666 0.34734356083861 0.38475913147671 0.42353138194053 0.46359738877080 0.504714692840

线性二步法近似值与精确值比较

2、四阶Runge-Kutta方法

function A=rungekutta() h=0.05; u(1)=0; for n=1:19 t(n)=h*(n-1); k1=2*t(n)-3*u(n);

k2=2*(t(n)+0.5*h)-3*(u(n)+0.5*h*k1); k3=2*(t(n)+0.5*h)-3*(u(n)+0.5*h*k2); k4=2*(t(n)+h)-3*(u(n)+h*k3); u(n+1)=u(n)+h/6*(k1+2*k2+2*k3+k4); end t(20)=1; A=[t' u']; plot(t',u','m');

legend('Runge-Kutta方法')

计算结果： ans =

0 0 0.05000000000000 0.00237968750000 0.10000000000000 0.00907095185669 0.15000000000000 0.01947322746655 0.20000000000000 0.03306960235350 0.25000000000000 0.04941516593756

0.30000000000000 0.06812697985904 0.35000000000000 0.08887544578091 0.40000000000000 0.11137687558199 0.45000000000000 0.135387095844 0.50000000000000 0.160695946779 0.55000000000000 0.18712253862403 0.60000000000000 0.21451118020302 0.65000000000000 0.24272786625619 0.70000000000000 0.27165725730431 0.75000000000000 0.30120007966637 0.80000000000000 0.33127088763203 0.85000000000000 0.36179613734408 0.90000000000000 0.392712527260 1.00000000000000 0.42396557981002

Runge-Kutta方法近似值与精确值比较：

精确解：

u=2/9*exp(-3*t)+2/3*(t-1/3) 0t1

精确值：

0 0 0.05000000000000 0.00237955031668 0.10000000000000 0.00907071570705 0.15000000000000 0.01947292258262 0.20000000000000 0.03306925246534 0.25000000000000 0.049414749800

0.30000000000000 0.06812659105347 0.35000000000000 0.08887505535803 0.40000000000000 0.111379153604 0.45000000000000 0.13538672458798 0.50000000000000 0.16069559114410 0.55000000000000 0.18712220191572 0.60000000000000 0.214510804924 0.65000000000000 0.24272757146367 0.70000000000000 0.27165698405622 0.75000000000000 0.30119982768041 0.80000000000000 0.33127065628654 0.85000000000000 0.36179592577803 0.90000000000000 0.392712336139 0.95000000000000 0.42396540463885 精确值曲线：