矩阵实验

有些实验可能不是很恰当，仅供参考

2.1

分别利用Gauss消元法和列主消去法编程求解方程组 Ax=b,其中

并求出矩阵A的LU分解及列主元的LU分解（求出L,U,P），并用LU分解的方法求A的逆矩阵及A的行列式

代码

 1function [x , L,U,Ani] = DelGauss(A,b)
 2% Gauss 消去法
 3% 输出解x,L,U,A的逆矩阵
 4Acopy=A;
 5n = length(b);
 6L = zeros(n,n);
 7for k = 1:n-1
 8    mul = A(k+1:n,k)/A(k,k);
 9    L(k+1:n,k) = mul;
10    A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - mul * A(k,k+1:n);
11    b(k+1:n) = b(k+1:n) - mul * b(k);
12    A(k+1:n,k) = zeros(n-k,1);
13end 
14x = zeros(n,1);
15x(n) = b(n)/ A(n,n);
16U = A;
17for k = n-1:-1:1
18    x(k) = (b(k) - A(k,k+1:n) * x(k+1:n))/A(k,k);
19end
20L(logical(eye(size(L))))=1;
21Uni=inv(U);
22Lni=inv(L);
23Ani=Uni*Lni;
24
25
26function [x,L,U,P] = ZGauss(A,b)
27% 列主元Guass 消去法
28% 输出解x,L,U,P
29n = length(b);
30L = zeros(n,n); 
31P = eye(n,n);
32for k=1:n-1
33    [ap,pos] = max(abs(A(k:n,k)));
34    pos = pos +k -1;
35    if pos > k
36        A([k,pos],:) = A([pos,k],:);
37        b([k pos],:) = b([pos,k],:);
38        P([k,pos],:) = P([pos,k],:);
39    end
40    mul = A(k+1:n,k) / A(k,k);
41    L(k+1:n,k) = mul;
42    A(k+1:n,k+1:n) = A(k+1:n,k+1:n) - mul * A(k ,k+1:n);
43    b(k+1:n) = b(k+1:n) - mul * b(k);
44    A(k+1:n,k) = zeros(n-k,1);
45end
46x = zeros(n,1);
47x(n) = b(n) / A(n,n);
48U = A;
49for k = n-1:-1:1
50    x(k) = (b(k) - A(k,k+1:n) * x(k+1:n)) / A(k,k);
51end
52L(logical(eye(size(L))))=1;

2.2

编制程序求解矩阵A的Cholesky分解，并用程序求解方程组Ax=b，其中

代码

 1function [X]=m_chol(A,b)
 2[N,N]=size(A);
 3X=zeros(N,1);
 4Y=zeros(N,1);
 5for i=1:N
 6    A(i,i)=sqrt(A(i,i)-A(i,1:i-1)*A(i,1:i-1)');
 7    if A(i,i)==0
 8           'A is singular. no unique solution'
 9        break
10    end
11    for j=i+1:N
12       A(j,i)=(A(j,i)-A(j,1:i-1)*A(i,1:i-1)')/A(i,i);
13    end
14end
15A
16b
17%前代法
18for j=1:N
19    Y(j)=(b(j)-A(j,1:j-1)*Y(1:j-1))/A(j,j);
20end
21Y
22%
23A=A'
24for k=N:-1:1
25    X(k)=(Y(k)-A(k,k+1:N)*X(k+1:N))/A(k,k);
26end

2.7

已知方程组

分别用Jacobi迭代法和Gauss-Seidel迭代法求解方程组，精确到小数点后6位，分别就n=10,20,30,50,100给出相应的计算结果

2.7.1 Jacobi迭代法

代码

 1function x=Jacobi(N)
 2%  x ---  方程组的解 
 3%  A --- 方程组的系数矩阵
 4A=diag(repmat([3], 1, N))+diag(repmat([-1], 1, N-1), 1)+diag(repmat([-1], 1, N-1), -1);
 5b=ones(N,1); % b --- 方程组的右端项
 6b(1)=2;
 7b(N)=2;
 8%disp(A);
 9%disp(b);
10x0 = zeros(length(b), 1) ;% x0 --- 解的初值 
11eps= 1.0e-6; % eps --- 精度要求，缺省为1e-6
12 
13M = 200;         % 最大迭代次数
14D=diag(diag(A)); % 求A的对角矩阵
15L=-tril(A,-1);   % 求A的下三角阵
16U=-triu(A,1);    % 求A的上三角阵
17B=D\(L+U); 
18f=D\b; 
19x=B*x0+f;  
20n=1; %迭代次数
21while norm(x-x0)>=eps 
22    x0=x;  
23    x=B*x0+f; 
24    n=n+1; 
25    if(n>=M)  
26        disp('Warning: ！'); 
27        return; 
28    end 
29end

2.7.2 Gauss-Seidel迭代法

代码

 1function [x]= Gauss_Seidel(N)
 2A=diag(repmat([3], 1, N))+diag(repmat([-1], 1, N-1), 1)+diag(repmat([-1], 1, N-1), -1);
 3b=ones(N,1); % b --- 方程组的右端项
 4b(1)=2;
 5b(N)=2;
 6    %y = zeros(1000,1);
 7    eps = 1.0e-6;
 8    D = diag(diag(A));
 9    L = -tril(A, -1);
10    U = -triu(A, 1);
11    B = (D - L)\U;
12    f = (D - L)\b;
13    count = 1;
14    x0 = zeros(N,1);
15    x = B*x0 + f;
16    tic;
17    while norm(x-x0) > eps
18        x0 = x;
19       % y(count) = norm(x-A\b);
20        x = B*x0 + f;
21        count = count + 1;
22
23        if count > 2000
24            disp('error:该矩阵不收敛');
25            return;
26        end
27    end
28    toc;
29    y = toc;
30    disp(count);
31end

2.8

用共轭梯度法求解上题中的方程组。

代码

 1function x=Gongetidu2(N)
 2%  x ---  方程组的解 
 3%  A --- 方程组的系数矩阵
 4A=diag(repmat([3], 1, N))+diag(repmat([-1], 1, N-1), 1)+diag(repmat([-1], 1, N-1), -1);
 5
 6b=ones(N,1); % b --- 方程组的右端项
 7b(1)=2;
 8b(N)=2;
 9%disp(A);
10%disp(b);
11x0 = zeros(length(b), 1) ;% x0 --- 解的初值 
12eps= 1.0e-6; % eps --- 精度要求，缺省为1e-6
13 
14M = 200;         % 最大迭代次数
15n=size(A,1);
16x=x0;
17r=b-A*x;
18d=r;
19for k=0:(n-1)
20    alpha=(r'*r)/(d'*A*d);
21    x=x+alpha*d;
22    r2=b-A*x;
23if ((norm(r2)<=eps)||(k==n-1))
24    x;
25    break;
26end
27beta=norm(r2)^2/norm(r)^2;
28d=r2+beta*d;
29r=r2;
30end
31end

3.4

用Newton迭代法求解方程x^{3}+x^{2}+x-3=0的根，初值选择x_{0}=-0.7,迭代7步并与真值x^{*}=1相比较，并列出数据表

从最后一列能得到什么结论？把最后的值与\frac{f^{''}(x^{*})}{2f^{'}(x^{*})}相比较能得到什么结论？

代码

 1function [x,ei]= Newton2()
 2syms x;
 3fun=inline('x^3+x^2+x-3');
 4dfun=inline('3*x^2+2*x+1');
 5x0=-0.7; 
 6N = 8;
 7ep = 1e-6; 
 8k=0;
 9while k<N 
10    x = x0-feval(fun,x0)/feval(dfun,x0) 
11    if abs(x-x0)<ep
12        break;
13    end
14    e = abs(x-1); %e=|x-x*| 1为真解x*
15    p = abs(e/(e^2-1));
16    k=k+1;
17    x0=x; 
18    ei = x-1
19end 
20if k==N
21    warning('已到达迭代次数上限！'); 
22end

3.5

取不同的初值用弦截法求方程x^{3}+2x^{2}+10x-100=0的实根，列表或者画图说明收敛速度。

代码

 1function [x,i]= xianjiefa(xx)%设置初值
 2t=[-10:0.1:10];
 3q1=t.^3+2*t.^2+10*t-100;
 4plot(t,q1);
 5%上 为此函数图像 可以看出函数解在3.4到3.5之间
 6%%%%%%%%%%%%%%%%%%%%%%%%
 7f = inline('x^3+2*x*x+10*x-100');
 8x = 0;
 9xxx = xx+1;
10i=0;
11while abs(xx-xxx) >= 1e-5
12	x = xx;
13	xx = xxx;
14	xxx = xx - f(xx) / (f(xx) - f(x)) * (xx - x);
15	i = i + 1;%i返回迭代次数
16end

3.6

设f(x)=54x^{6}+15x^{5}-102x^{4}-69x^{3}+35x^{2}+16x-4。在区间[-2,2]上画出函数，(1)使用Newton迭代法找出该区间上的5个根，并计算e_{i+1}/e^{2}_{i}和e_{i+1}/e_{i}，由此判断哪个根是1阶收敛，哪个根是2阶收敛？(2)使用分割线法计算这5个根，并判断哪个根是线性收敛，哪个是超线性收敛？

3.6.1 Newton迭代法

代码

 1t=[-2:0.01:2];
 2q1=54*t.^6+45*t.^5-102*t.^4-69*t.^3+35*t.^2+16*t-4;
 3plot(t,q1);%画出该函数图像
 4%%%%%%%%%%%%%%%%%
 5function [k,xk,yk,piancha,i]=newtonqx(x0,tol,gmax)
 6global fnq dfnq
 7x(1)=x0;
 8for i=1:gmax
 9    x(i+1)=x(i)-fnq(x(i))/(dfnq(x(i)+eps)); 
10    piancha = abs(x(i+1)-x(i));   
11    i=i+1;
12    xk=x(i)  
13    yk=fnq(x(i));   
14    if( piancha<tol ) 
15        k=i-1;
16        xk=x(i);
17        yk=fnq(x(i));   
18        break;
19    end
20end
21if i>gmax  
22    disp('超过最大迭代次数')
23    k=i-1;
24    xk=x(i);
25    yk=fnq(x(i));
26    %[i-1 xk yk piancha];
27    return;
28end

3.6.2 割线法

代码

 1function [k,xk,yk,piancha,i]=gexian(x01,x02,tol,gmax)
 2global fnq dfnq
 3x(1)=x01;
 4x(2)=x02;
 5for i=2:gmax
 6    %割线法迭代公式:  x(n+1)=x(n)- f(x(n))*( x(n)-x(n-1))/(f(x(n))-f(x(n-1) )
 7    % 即用x(n),x(n-1)上的差商替代导数f'(xn)
 8    u(i)=fnq(x(i))*(x(i)-x(i-1));
 9    v(i)=fnq(x(i))-fnq(x(i-1));
10    x(i+1)=x(i)-u(i)/(v(i));
11  
12    piancha=abs(x(i+1)-x(i));  %计算迭代精度
13
14    i=i+1;  
15    xk=x(i);
16    yk=fnq(x(i));
17  
18    if(piancha<tol)
19        k=i-2;
20        xk=x(i);
21        yk=fnq(x(i));
22        return  
23    end
24end

4.3

令f(x)= e^{|x|},x \in [-1,1]，分别用等距节点和Chebyshev的零点去插值f(x)，等距节点包括左右两个端点，分别取n= 5,10,15,20，画出插值函数以及原函数的图并比较，观察有没有Runge现象发生。

代码

lagrange

 1function y=lagrange(x0,y0,x)
 2ii=1:length(x0); 
 3y=zeros(size(x));
 4for i=ii
 5    ij=find(ii~=i);   
 6    y1=1;
 7    for j=1:length(ij)
 8        y1=y1.*(x-x0(ij(j)));  
 9    end
10    y=y+y1*y0(i)/prod(x0(i)-x0(ij));
11end

chebyshev

1function x=chebyshev(n)
2x = zeros(1,n+1);
3for k = 1:n+1
4    x(k) = cos(((2*k-1)*pi)/(2*(n+1)));
5end

omega

1function err = omega(x,x0)
2err = zeros(1,length(x));
3for i = 1:length(x)
4    err(i) = prod(x(i)*ones(1,length(x0))-x0);
5end
6err = err*exp(1)/factorial(length(x0));

测试程序

 1function main()
 2    x = -1:0.01:1;
 3    y = exp(1).^abs(x);
 4    %均匀间距插值，使用等间距节点作为插值节点
 5    subplot(3,2,1);
 6    x01 = -1:0.4:1; %n=5
 7    y01 = exp(1).^abs(x01);
 8    f01 = lagrange(x01,y01,x);
 9    err01 = exp(1)*omega(x,x01);
10    plot(x,f01,'-c')
11    title('均匀间距插值');
12    hold on
13    x11 = -1:0.2:1; %n=10
14    y11 = exp(1).^abs(x11);
15    f11 = lagrange(x11,y11,x);
16    err11 = exp(1)*omega(x,x11);
17    plot(x,f11,'-r')
18    hold on
19    x02 = -1:0.133:1; %n=15
20    y02 = exp(1).^abs(x02);
21    f02 = lagrange(x02,y02,x);
22    err02 = exp(1)*omega(x,x02);
23    plot(x,f02,'-g')
24    hold on
25    x12 = -1:0.1:1; %n=20
26    y12 = exp(1).^abs(x12);
27    f12 = lagrange(x12,y12,x);
28    err12 = exp(1)*omega(x,x12);
29    plot(x,f12,'-b')
30    hold on
31    plot(x,y,'-.k')
32    legend('n=5','n=10','n=15','n=20','原函数')
33
34    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%5
35    %切比雪夫插值，使用切比雪夫节点作为插值节点
36    subplot(3,2,2);
37    x23 = chebyshev(5);%n=5
38    y23 = exp(1).^abs(x23);
39    f23 = lagrange(x23,y23,x);
40    err23 = exp(1)*omega(x,x23);
41    plot(x,f23,'-c');
42    hold on
43    x21 = chebyshev(10);%n=10
44    y21 = exp(1).^abs(x21);
45    f21 = lagrange(x21,y21,x);
46    err21 = exp(1)*omega(x,x21);
47    plot(x,f21,'-r');
48    hold on
49    x24 = chebyshev(15);%n=15
50    y24 = exp(1).^abs(x24);
51    f24 = lagrange(x24,y24,x);
52    err24 = exp(1)*omega(x,x24);
53    plot(x,f24,'-g');
54    x22 = chebyshev(20);%n=20
55    y22 = exp(1).^abs(x22);
56    f22 = lagrange(x22,y22,x);
57    err22 = exp(1)*omega(x,x22);
58    plot(x,f22,'-b');
59    plot(x,y,'-.k')
60    title('切比雪夫插值');
61    legend('n=5','n=10','n=15','n=20','原函数')
62    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
63    %绘制插值余项的最大值
64    subplot(3,2,3);
65    plot(x,err01,'-r',x,err23,'-b');
66    title('当为5次插值余项最值');
67    legend('等距节点','切比雪夫节点');
68    subplot(3,2,4);
69    plot(x,err11,'-r',x,err21,'-b');
70    title('当为10次插值余项最值');
71    legend('等距节点','切比雪夫节点');
72    subplot(3,2,5);
73    plot(x,err02,'-r',x,err24,'-b');
74    title('当为15次插值余项最值');
75    legend('等距节点','切比雪夫节点');
76    subplot(3,2,6);
77    plot(x,err12,'-r',x,err22,'-b');
78    title('当为20次插值多项式时的插值余项最值');
79    legend('等距节点','切比雪夫节点');

4.4

(1)给定数据点(x_{i},x^{2}_{i}),x_{i}-0,\frac{1}{n},\frac{2}{n},\cdots ,1，当n=5,10,15,20,25,30时分别用直线拟合这组数据点并注意观察当点数逐渐增加时直线的表达式的变化。(2)计算函数f\left( c_{1},c_{2}\right) =\int^{1}_{0} \left( x^{2}-c_{1}-c_{2}x\right)^{2} dx的最小值，并解释与(1)的关系

代码

 1function  Least_square(N)
 2    x =(0:1/N:1);
 3    y=x.*x;
 4    P=polyfit(x,y,1);
 5    y001=polyval(P,x);  
 6    plot(x,y001,x,y,'r*'); 
 7    title(['n=',num2str(N)]);
 8%%%%%%%%%%%%
 9Least_square(5)
10Least_square(10)
11Least_square(15)
12Least_square(20)
13Least_square(25)
14Least_square(30)

4.6

代码

 1function [f,RMSe] = Least_square(order)%参数为阶数
 2    x=[1994 1995 1996 1997 1998 1999 2000 2001 2002 2003]; 
 3    y=[67.052 68.008 69.803 72.024 73.400 72.063 74.669 74.487 74.065 76.777];  
 4    P=polyfit(x,y,order);
 5    x1=1994:0.1:2003;  
 6    y1=polyval(P,x1);  
 7    plot(x1,y1,x,y,'r*'); 
 8    xlabel('年份');
 9    ylabel('百万桶/天（\times10^6）');
10  
11    x2=1994:1:2003; 
12    y2=polyval(P,x2);
13    RMSe=sqrt(sum((y2-y).^2)/10);
14    f=vpa(poly2sym(P));
15%%%%%%%%
16%执行代码
17 [f,RMSe] = Least_square(3)
18 [f,RMSe] = Least_square(2)
19 [f,RMSe] = Least_square(1)

4.7

编程计算三次样条S，满足S(0)=1,S(1)=3,S(2)=3,S(3)=4,S(4)=2,其中边界条件S^{''}(0)=S^{''}(4)=0。

代码

 1function cubic_spline()
 2    % x=input('请按照格式[x1,x2,x3...]格式输入y=f(x)函数已知点的横坐标xi='); 
 3    % y=input('请按照格式[y1,y2,y3...]格式输入y=f(x)函数已知点对应的纵坐标yi=');
 4    x = [0,1,2,3,4];
 5    y = [1,3,3,4,2];
 6    n=size(x,2); 
 7    for k=2:n 
 8    h(k-1)=x(k)-x(k-1); 
 9    end
10    for k=1:(n-2) %计算μ和λ 
11    mu(k)=h(k)/(h(k)+h(k+1)); 
12    lambda(k)=1-mu(k); 
13    end
14    for k=1:(n-2) 
15    g(k)=3*(lambda(k)*(y(k+1)-y(k))/h(k)+mu(k)*(y(k+2)-y(k+1))/h(k+1)); %计算g(1)到g(n-2) 
16    end   
17    %     M(1)=input('请输入f(a)的二阶导数值：'); 
18    %     M(n)=input('请输入f(b)的二阶导数值：');  
19        M(1)=0; 
20        M(n)=0;  
21        A=zeros(n,n);     %构造追赶法所需的A和b 
22    for k=2:(n-1)   
23            A(k,k)=2; 
24            A(k,k+1)=mu(k-1); 
25            A(k,k-1)=lambda(k-1); 
26    end 
27
28    A(1,1)=2; 
29        A(1,2)=1; 
30        A(end,end)=2; 
31        A(end,end-1)=1; 
32        A; 
33        b=zeros(n,1);
34
35    for k=2:(n-1) 
36            b(k,1)=g(k-1); 
37    end
38
39        b(1,1)=3*((y(2)-y(1))/h(1)-2*h(1)*M(1)); 
40        b(n,1)=3*((y(n)-y(n-1))/h(n-1)+2*h(n-1)*M(n)); 
41        b; 
42        b=b'; 
43        m=chase(A,b);  %利用追赶法求解成功，这里的参数b形式应为行向量而非列向量  
44
45    for k=1:(n-1) 
46      clear S1 
47      syms X 
48      S1=(1-2*(X-x(k))/(-h(k)))*((X-x(k+1))/(h(k)))^2*y(k)+(X-x(k))*((X-x(k+1))/(h(k)))^2*m(k)+(1-2*(X-x(k+1))/(h(k)))*((X-x(k))/(h(k)))^2*y(k+1)+(X-x(k+1))*((X-x(k))/(h(k)))^2*m(k+1);
49      fprintf('当%d=<X=<%d时\n',x(k),x(k+1)); 
50      S=expand(S1)
51    end
52%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
53%追赶法
54function [x] = chase(A,d)
55    if ~isequal(tril(A,-1)-diag(diag(tril(A,-1),-1),-1),zeros(size(A)))
56    error('输入矩阵不是一个三角矩阵');
57    end
58    if ~isequal(triu(A,1)-diag(diag(triu(A,1),1),1),zeros(size(A)))
59    error('输入矩阵不是一个三角矩阵');
60    end
61    a=[0;diag(tril(A,-1),-1)];%下对角线
62    b=diag(A);%中对角线
63    c=[diag(triu(A,1),1);0];%上对角线
64    l=zeros(size(a,1),1);%求L
65    u=zeros(size(b,1),1);%求U
66    n=size(b,1);%矩阵的维度
67    x=zeros(n,1);
68    y=zeros(n,1);
69    u(1)=b(1);
70    for i=2:n
71        l(i)=a(i)/u(i-1);
72        u(i)=b(i)-l(i)*c(i-1);
73    end
74    y(1)=d(1);
75    for i=2:n
76        y(i)=d(i)-l(i)*y(i-1);
77    end
78    x(n)=y(n)/u(n);
79    for i=n-1:-1:1
80        x(i)=(y(i)-c(i)*x(i+1))/u(i);
81    end
82end
83%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
84%生成图像
85x=0:0.1:1;
86x1=1:0.1:2;
87x2=2:0.1:3;
88x3=3:0.1:4;
89hold on
90plot(x, (149*x)/56 - (37*x.^3)/56 + 1)
91plot(x1, (73*x1.^3)/56 - (165*x1.^2)/28 + (479*x1)/56 - 27/28)
92plot(x2, - (87*x2.^3)/56 + (45*x2.^2)/4 - (1441*x2)/56 + 613/28)
93plot(x3, (51*x3.^3)/56 - (153*x3.^2)/14 + (2285*x3)/56 - 625/14)

5.1

已知f(x)=x^{2}\sin x，分别用复化梯度形公式和复化Simpson公式计算积分\int^{pi }_{0} f(x)dx，区间分别为20,40,80,200个小区间，并计算其精确值，比较计算精度情况。

5.1.1 精确值

5.1.2 复化梯形公式

代码

 1function [y_x] = CalcuFunctionValue(x)
 2if x == 0
 3    y_x = 1;
 4else
 5    y_x = sin(x)*x*x;
 6end
 7end
 8
 9function [res] = ComplexTrap( x_low, x_up,n)
10step_length = (x_up - x_low)/n;
11res = 0;
12for i = 1:n-1
13    res = res + CalcuFunctionValue(x_low+step_length*(i-1))+CalcuFunctionValue(x_low+step_length*i);
14end 
15res = res * step_length / 2;
16end

5.1.3 复化Simpson公式

代码

 1function [res] = ComplexSimpson(x_low, x_up,n)
 2% simpson求积公式
 3if mod(n,2) == 0
 4    step_length = (x_up - x_low)/n;
 5    res = 0;
 6    for i = 1:2:n-1
 7        res = res + CalcuFunctionValue(x_low+step_length*(i-1))...
 8                 +4*CalcuFunctionValue(x_low+step_length*i)...
 9                 +CalcuFunctionValue(x_low+step_length*(i+1));
10    end 
11    res = res*step_length/6*2;
12else
13    print('等分区间数错误！');
14end 
15end

6.1

已知常微分方程

分别用Euler法，改进的Euler法，Runge-Kutta法去求解该方程，步长选为0.1,0.05,0.01.画图观察求解效果。

代码

Euler

 1function [y_x] = CalcuFunctionValue(x,y) %所求的函数
 2if x == 0
 3    y_x = 1;
 4else
 5    y_x = 2/x*y+x*x*exp(x);
 6end
 7end
 8
 9function [x, y] = forwardEuler(x0, y0, x1, h)
10n = floor((x1 - x0)/h);
11x = zeros(n + 1, 1);
12y = zeros(n + 1, 1);
13x(1) = x0;
14y(1) = y0;
15for i = 1 : n
16    x(i + 1) = x(i) + h;
17    y(i + 1) = y(i) + h * CalcuFunctionValue(x(i),y(i));
18end

改进的Euler

 1function [x, y] = improvedEuler(x0, x1, y0, h)
 2n = floor((x1 - x0)/h);
 3x = zeros(n + 1, 1);
 4y = zeros(n + 1, 1);
 5x(1, 1) = x0;
 6y(1, 1) = y0;
 7for i = 2 : n+1
 8    x(i, 1) = x(i - 1, 1) + h;
 9    y(i, 1) = y(i - 1, 1) + h *CalcuFunctionValue(x(i - 1, 1),y(i - 1, 1));
10    y2 = y(i - 1, 1) + h  *  CalcuFunctionValue(x(i, 1),y(i, 1));
11    y(i, 1) = (y(i, 1)+y2)/2
12end

Runge-Kutta

 1function [x, y] = runge(x0, x1, y0, h)
 2n = (x1 - x0) / h;
 3x = zeros(n + 1);
 4y = zeros(n + 1);
 5x(1) = x0;
 6y(1) = y0;
 7for i = 1:n
 8    x(i + 1) = x(i) + h;
 9    k1 = CalcuFunctionValue(x(i), y(i));
10    k2 = CalcuFunctionValue(x(i) + 0.5*h, y(i) + k1*h/2);
11    k3 = CalcuFunctionValue(x(i) + 0.5*h, y(i) + k2*h/2);
12    k4 = CalcuFunctionValue(x(i)+ h, y(i) + k3*h);
13    y(i + 1) = y(i) + h*(k1 + 2*k2 + 2*k3 + k4)/6;
14end
15end

P110例4 P243例1

• 这个题理解的不太好，不知道做的对不对

代码

精细算法

 1function x = jingxi(A,x0,h)
 2    n=length(x0);
 3    I=eye(n);
 4    x(1:n,1)=x0;
 5    N=20;
 6    dt=h/2^N;
 7    At=A*dt;
 8    BigT=At*(I+At*(I+At/3*(I+At/4))/2);
 9    for k=1:N
10        BigT=2*BigT+BigT^2;
11    end
12        BigT = I + BigT;
13    for k=1:n
14        x(:,k+1)=BigT*x(:,k);
15    end

四阶Runge-Kutta算法

 1function F= f(t,Y)
 2x=Y(1);
 3y=Y(2);
 4f1=-2000*x+999.75*y+1000.25;
 5f2=x-y;
 6F=[f1;f2];
 7end
 8%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 9function Rk4_test(Delta,x_zero,y_zero)
10   % Delta =0.001;
11    t=0:Delta:8;
12    n=length(t);
13    Y(:,1)=[x_zero;y_zero];
14    for k=1:n-1
15        z1 = f(t(k),Y(:,k));
16        z2 = f(t(k)+Delta/2,Y(:,k)+z1*Delta/2);
17        z3 = f(t(k)+Delta/2,Y(:,k)+z2*Delta/2);
18        z4 = f(t(k)+Delta,Y(:,k)+z2*Delta);
19        Y(:,k+1) = Y(:,k)+Delta*(z1+2*z2+2*z3+z4)/6;
20    end
21    x=Y(1,:);
22    y=Y(2,:);
23
24    figure
25    subplot(1,2,1)
26    plot(t,x,'b')
27    xlabel('t')
28    ylabel('x')
29    subplot(1,2,2)
30    plot(t,y,'r')
31    xlabel('t')
32    ylabel('y')
33%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
34%test
35Rk4_test(0.001,0,0)

总结

写了持续一个多周的数值实验，本来想用python写，但没用过MATLAB，就像用用，原来这个软件这么强大，速成了MATLAB，收获很大，可能有些数值实验题目做的不是很对，但终归是自己借助网络和自己敲敲改改出来的。

截图

实验报告截图链接

    

公众号'艾恩凝' 个人微信

    吾心信其可行，
          则移山填海之难，
                  终有成功之日！
                                  ——孙文