科学计算实践上机实验报告

matlab 计算

第一次上机作业

题目

用积分方程的离散解并比较它与真解之间的误差
用 $n=4$ 的复化公式求积分方程

$$\frac{4x^3+5x^2-2x+5}{8(x+1)^2}+\int \limits^1_0 (\frac{1}{1+t}-x)y(t)dt=y(x)$$ 的近似解，并在区间 $[0,1]$ 的一些点上与真解 $y(x)=\frac{1}{(1+x)^2}$ 进行比较

分析

• 由题目分析可知，未知函数为 $y(x)$ ,利用复化辛普森求积公式 $$\int \limits_a^bf(x) = \frac{h}{6}[f(a)+2\sum\limits^{n-1}_{i=1}f(x_i)+4\sum\limits_{i=1}^{n}f(\frac{x_{i-1+x_i}}{2})+f(b)$$ 可以将方程的中的积分形式变成离散的格式
• 当取 $n$ 个插值点时，此时对积分区间$[0,1]$ 应取的插值解点为 $2*n+1$ 个。所以 $n=4$ 插值节点个数为9个
• 对于离散的方程中，因为积分的结果中仍然包含有未知数$x$ ， 因此不能直接获得未知方程$y(x)$ 在离散点处的值
• 此时将离散方程中 $x$ 取值为辛普森积分中的插值点，这样就获得了离散的方程组。
• 对方程组进行求解便可获得未知方程 $y(x)$ 在插值节点处的值
• 使用三次样条插值对$y(x)$ 进行插值并画出图像与函数的真解进行对比

代码

主函数
代码如下：

%% 
n = 6;
h = 1/n;
x =0:h/2:1;
a=ones(1,(2*n+1));
fre = ones(1,(2*n+1));
%% 这段用来定义符号变量 获得 2*n 个符号变量并存入y 以向量的方式提取
cha = '[ ';
for i=1:(2*n+1)
     cha = strcat(cha,' y',num2str(i)) ;   
end
cha = strcat(cha,']');
y = sym(cha);
%%
%%定义第一个函数，    
f1 = inline('(4*(x^3)+5*(x^2)-2*x+5)/(8*((x+1)^2))','x'); 
%% 算出方程组每个方程的表达式
solveequa = {};
for i=1:2*n+1;
    a(i) = f1(x(i));
    rs= a(i)+sym(fuhuaxin(n,x(i),y))-y(i);
    resu = char(rs);
    %resul = char(y(1));
    res = strcat(resu,'= 0');
    solveequa{i} = (res );
end
%%求解方程并将方程的解存为向量的形式
%y= solve(solveequa{1},solveequa{2},solveequa{3},solveequa{4},solveequa{5},solveequa{6},solveequa{7},solveequa{8},solveequa{9});
y = solve(solveequa{:});
y = struct2cell(y);
yang = 1:2*n+1 ;
for i=1:2*n+1;
    yang(i) = vpa(y{i});
end
%% 三次样条差值
newx = 0:0.1:1;
simy = interp1(x,yang,newx,'spline');
%% 绘出差值的图像
plot(newx,simy,'-*');
hold on
f3 = inline('1/((x+1)^2)','x');
%% 绘出真实的函数值的图像
x = 0:0.01:1;
yang= x;
for i = 1:101;
    yang(i) = f3(x(i));
end
plot(x,yang,'r');

• 其中子函数fuhuaxin.m
  代码如下：

function[res] = fuhuaxin(n,x,y)
f2 = inline('(1/(t+1)-x)*y','x','t','y');
h = 1/n;
t = 0:h/2:1;
res = f2(x,t(1),y(1));
for i=1:(n-1)
    res = res + 2*f2(x,t(2*i+1),y(2*i+1));
end
for i =1:n
    res = res +4* f2(x,t(2*i),y(2*i));
end
res = res + f2(x,t(2*n+1),y(2*n+1));
res = (h/6)*res;

结果

• 各个点处的数值为：
• 函数图像
  
  • 其中蓝色表示利用复化辛普森方程和三次样条插值求出来的未知函数 $y(x)$ 图像
  • 其中红色表示函数真解的图像
  • 根据图像分析可以看出来，两个函数相差很小，因此可以认为所求的的插值解是有效的解

第二次上机作业

第一题

题目

分别用共轭梯度算法和预条件共轭梯度算法求解线性方程组 $Ax=b$ ，其中

$$A=\begin{pmatrix}4&1\\1&4^2&1\\&1&4^3&1\\&&1&4^4&1\\&&&1&4^5&1\\&&&&1&4^6&1\\&&&&&1&4^7&1\\&&&&&&1&4^8&1\\&&&&&&&1&4^9&1\\&&&&&&&&1&4^{10}&1\\\end{pmatrix},b=\begin{pmatrix}5\\8\\66\\258\\1026\\4098\\16386\\65538\\262146\\1048577\end{pmatrix}$$

分析

• 因为矩阵形式为对成正定且对角占优，因此可以使用共轭梯度法求解线性方程组，需要选定初始值，下降方向和下降步长。
• 由共轭梯度，给定初始向量 $x_0$ ，第一步使用负梯度方向为下降方向,即 $p_0=r_0$ ,于是有 $$\alpha_0=\frac{r_0^Tr_0}{p_0^TAp_0},\quad x_1=x_0+\alpha_0p_0,\quad r_1=b-Ax_1$$
• 对于以后各部，如 $k+1$ 步$(k\geq1)$ ,下降方向就不再取 $r_k$ ,而是 $$\alpha_k=\frac{r_k^Tp_k}{p_k^TAp_k}\\\beta_K=-\frac{r_{k+1}^TAp_k}{p_k^TAp_k}$$
• 当 $$x_{k+1}-x_k\leq \varepsilon$$ 时，使迭代停止。这时的$x_k$就是所求的方程组的解

代码

主函数代码如下 GG_main：

function GG_main()
G = zeros(10)
for t = 1:10
    G(t,t) = 4^t;
    if t>=2
        G(t,t-1) = 1;
    end
    if t<=9
        G(t,t+1) = 1;
    end
end
G
b =[5 18 66 258 1026 4098 16386 65538 262146 1048577]
x0 = zeros(1,10)
max_iter = 1000;
fprintf('\n')
fprintf('Conjugate Gradient Methord:\n');
fprintf('=============================');
[y,iter] = cg(G, b, x0, max_iter );
fprintf('\n')
fprintf(' Iterative number: \n      %d \n',iter);
fprintf(' Solution:\n');
fprintf('   %10.6f',y);
fprintf('\n\n==========================\n')

cg 函数代码：

function [x iter] = cg(G, b, x0, max_iter)
x = x0';
b=b'
tolerance = 1.0e-6;
fprintf('\n x0=');
fprintf('   %10.6f',x0);
%save sb.mat G b x
r = G*x-b;
d=-r;
for k=1:max_iter
    if norm(r,2) <= tolerance
        iter = k-1;
        fprintf('\n Algorithm finds a solution!');
        return
    end
    alpha = (r'*r)/(d'*G*d);
    xx = x+alpha*d;
    rr = r+alpha*G*d;
    beta = (rr'*rr)/(r'*r);
    d = -rr+beta*d;
    x = xx;
    r = rr;
    fprintf('\n  x%d = ',k);
    fprintf('    %10.6f',x);
end
iter = max_iter;
return 

结果

• 方程原型：
• 迭代过程：
• 观察迭代结果可知：迭代所获得的解与方程组的真解几乎没有误差；因此可以认为共轭梯度法迭代方法是有效的

第二题

题目

用预条件共轭梯度算法求解线性方程组 $Ax=b$ , 其中

$$A=\begin{pmatrix}B &-I\\-I&B&-I\\&\ddots&\ddots&\ddots\\&&-I&B&-I\\&&&-I&B\end{pmatrix},\quad b=\begin{pmatrix}c\\d\\\vdots\\d\\c\end{pmatrix}$$ 上式中 $A$ 为 $N^2\times N^2$ 维矩阵，$b$ 为 $N^2$ 维向量， $B$ 为 $N\times N$ 矩阵，$I$ 为 $N\times N$ 单位矩阵，$c,d$ 为 $N$ 维向量，且 $$B=\begin{pmatrix}4 &-1\\-1&4&-1\\&\ddots&\ddots&\ddots\\&&-1&4&-1\\&&&-1&4\end{pmatrix},\quad c=\begin{pmatrix}2\\1\\\vdots\\1\\2\end{pmatrix},\quad b=\begin{pmatrix}1\\0\\\vdots\\0\\1\end{pmatrix},\quad N=400$$

分析

• 分析题目可知，线性方程组的规模超过了计算机的存储范围，因此不能在运算过程中直接带入矩阵或是向量，因而需要使用 crs 压缩行存储的方法存储矩阵
• 在运算过程中也是使用矩阵中非零元素，这样可以减小运算规模便于求解。
• 使用预优共轭梯度法，因此需要选定预优矩阵，因为矩阵的规模较大，而且对角线元素相同，都为 4 ,所以选定与有矩阵为对角矩阵，元素为 4
• 因为矩阵形式为对成正定且对角占优，因此可以使用共轭梯度法求解线性方程组，需要选定初始值，下降方向和下降步长。
• 由共轭梯度，给定初始向量 $x_0$ ，第一步使用负梯度方向为下降方向,即 $p_0=r_0$ ,于是有 $$\alpha_0=\frac{r_0^Tr_0}{p_0^TAp_0},\quad x_1=x_0+\alpha_0p_0,\quad r_1=b-Ax_1$$
• 对于以后各部，如 $k+1$ 步$(k\geq1)$ ,下降方向就不再取 $r_k$ ,而是 $$\alpha_k=\frac{r_k^Tp_k}{p_k^TAp_k}\\\beta_K=-\frac{r_{k+1}^TAp_k}{p_k^TAp_k}$$
• 当 $$x_{k+1}-x_k\leq \varepsilon$$ 时，使迭代停止。这时的$x_k$就是所求的方程组的解
• 求出方程组的解后，因为方程组的规模过大，所以解输出到文件当中查看。

代码

代码主要有四个函数 main cg getmatrix getvalue 其中
1. main 为运算的入口；

function [] = main()
%% 请输入N
N = input('Please input N:');
%% 通过getmatrix 函数获得矩阵
[b, value, cow, row] = getmatrix(N);
%% 获得向量
max_it = 1000;
% for m=1:N*N
%     for n=1:N*N
%         vae = getvalue(m,n,value,cow,row);
%         fprintf('%g\t',vae);
%     end
%     fprintf('\n')
% end
% disp(b);
M = 1/4;
x0 = zeros(N*N,1);
[x, iter] = cg(b, x0,M, max_it, value, cow, row);
xlswrite('test1.xlsx',x,1);
iter
end

1. getmatrix 根据 $N$ 的大小获得相应的矩阵 $A$ 的 crs 压缩行存储格式

function [b, value, cow, row]=getmatrix(N)
% 这个函数的作用是产生矩阵的crs 存储 返回值是 值向量value， 列向量cow 以及行向量row。
%% 对向量开始初始化
valuelen = (N*3-2)*N + N*(N-1)*2;
value = ones(valuelen,1);
cow = ones(N*N,1);
row = ones(valuelen,1);
tempadr = 1;%利用 tempadr 记录 value 数值变化的位置;
cowadr = 2;%利用 cowadr 记录 cow 列号变化的位置；
rowadr = 1 ;%利用 row 记录行变化的位置
%% 开始形成b
c = ones(N,1);
c(N) = 2;
c(1) = 2; 
d = zeros(N,1);
d(N) = 1;
d(1) = 1;
b = zeros(N*N,1);
for t = 0:N-1
    temp = t*N+1;
    endp = (t+1)*N;
    if t>0 && t<N-1
        b(temp:endp) = d;
    else
        b(temp:endp) = c;
    end
end
    %% 开始计算
for t=1:N
    if t==1
        for j=1:N
            if j==1
            value(tempadr:tempadr+2) = [4 -1 -1];
            tempadr = tempadr + 3;
            cow(cowadr) = tempadr;
            cowadr = cowadr +1;
            row(rowadr:rowadr+2) =[1 2 N+1];
            rowadr = rowadr + 3;
            else
                if j<N                    
                    value(tempadr:tempadr+3) = [-1 4 -1 -1];
                    tempadr =tempadr + 4;
                    cow(cowadr) = tempadr;
                    cowadr = cowadr +1;
                    row(rowadr:rowadr+3) =[j-1 j j+1 j+N];
                    rowadr = rowadr + 4;
                else
                    value(tempadr:tempadr+2) = [-1 4 -1];
                    tempadr = tempadr + 3;
                    cow(cowadr) = tempadr;
                    cowadr = cowadr +1;
                    row(rowadr:rowadr+2) =[j-1 j j+N];
                    rowadr = rowadr + 3;
                end
            end
        end   
    else           
        if t<N
            for j=1:N
                if j==1
                value(tempadr:tempadr+3) = [-1 4 -1 -1];
                tempadr = tempadr + 4;
                cow(cowadr) = tempadr;
                cowadr = cowadr +1;
                row(rowadr:rowadr+3) =[(t-2)*N+j (t-1)*N+j (t-1)*N+j+1  t*N+j];
                rowadr = rowadr + 4;
                else
                    if j<N                    
                        value(tempadr:tempadr+4) = [-1 -1 4 -1 -1];
                        tempadr = tempadr + 5;
                        cow(cowadr) = tempadr;
                        cowadr = cowadr +1;
                        row(rowadr:rowadr+4) =[(t-2)*N+j (t-1)*N+j-1 (t-1)*N+j (t-1)*N+j+1 t*N+j];
                        rowadr = rowadr + 5;
                    else
                        value(tempadr:tempadr+3) = [-1 -1 4 -1];
                        tempadr = tempadr + 4;
                        cow(cowadr) = tempadr;
                        cowadr = cowadr +1;
                        row(rowadr:rowadr+3) =[(t-2)*N+j (t-1)*N+j-1 (t-1)*N+j  t*N+j];
                        rowadr = rowadr + 4;
                    end
                end
            end
        else
            for j=1:N
                if j==1
                value(tempadr:tempadr+2) = [-1 4 -1];
                tempadr = tempadr + 3;
                cow(cowadr) = tempadr;
                cowadr = cowadr +1;
                row(rowadr:rowadr+2) =[(t-2)*N+j (t-1)*N+j (t-1)*N+j+1 ];
                rowadr = rowadr + 3;
                else
                    if j<N                    
                        value(tempadr:tempadr+3) = [-1 -1 4 -1 ];
                        tempadr = tempadr + 4;
                        cow(cowadr) = tempadr;
                        cowadr = cowadr +1;
                        row(rowadr:rowadr+3) =[(t-2)*N+j (t-1)*N+j-1 (t-1)*N+j (t-1)*N+j+1];
                        rowadr = rowadr + 4;
                    else
                        value(tempadr:tempadr+2) = [-1 -1 4 ];
                        tempadr = tempadr + 3;
                        cow(cowadr) = tempadr;
                        cowadr = cowadr +1;
                        row(rowadr:rowadr+3) =[(t-2)*N+j (t-1)*N+j-1 (t-1)*N+j (t-1)*N+j+1];
                        rowadr = rowadr + 4;
                    end
                end
            end
        end            
    end
end

1. getvalue 根据压缩行存储格式获得矩阵相应行和列的数值（主要用来测试生成的格式是否正确）

function[val] = getvalue(m,n,value,cow,row)
% 这个函数通过 crs 存储的 value , cow ,row 计算出矩阵第m行 n列的元素的值并返回
%%
[lenth,~] = size(cow);
if m>= lenth
     fprintf('Larger than the matrix\n')
else
    val = 0;    
    for t=cow(m):(cow(m+1)-1)  
        if n==row(t)        
            val = value(t);
            return
        end
    end
end
end

1. cg 进行预优共轭梯度迭代的函数

function [x, iter] = cg(b, x0,M, max_it, value, cow, row)
x=x0;
iter = 0;
tolerance = 1.0e-6;
N = length(b);
%% 现在计算 A*x 的值
temp= zeros(N,1);
for t = 1:N
    for j = cow(t):(cow(t+1)-1) 
        temp(t) = temp(t) + value(j)*x(row(j));   
    end
end
r = b - temp;
%% 开始循环
for k=1:max_it
    if norm(r,2) <= tolerance*norm(b)
        iter = k-1;
        fprintf('\n Algorithm finds a solution!');
        return
    end
    z = M*r;
    iter = iter + 1;
    if iter==1
        p = z;
        roo = r'*z;
    else
        roos = roo;
        roo = r'*z;
        beta = roo*(1/roos);
        p = z + beta*p;
    end
    temp= zeros(N,1);
    for t = 1:N
        for j = cow(t):(cow(t+1)-1) 
            temp(t) = temp(t) + value(j)*p(row(j));   
        end
    end
    w = temp;
    alpha = roo*(1/(p'*w));
    x = x + alpha*p;
    r = r - alpha * w;
end
end

结果

• 取 $N=3$ 将矩阵输出查看 crs压缩行存储的格式是否正确：
  可以看出数据的格式正确，最终方程的解为全部为 1 的列向量，长度为 $N\times N$ .
• 因此预优共轭梯度法结合 crs 压缩行存储格式迭代程序是有效的。

第三次上机作业

第一题

题目

用差分法求解如下二阶常微分方程边值问题 $(N=400)$

$$\left\{ \begin{array}{rcl} -\frac{d^2u}{dx^2}=x^2&&0<x<1,\\ u(0)=1,&u'(1)+2u(1)=1. \end{array} \right.$$

分析

• 差分法求解二阶常微分方程，根据题目条件可知，边值条件为第三边值条件
• 不妨记 $f(x) = -x^2$，我们根据差分方程 $$y_{i-1}-(2+h^2q_i)y_i+y_{i+1}=h^2f_i,\quad i=1,2,\cdots,n-1$$
• 获得的差分方程组为 $$\left\{\begin{array}{rcl}-(3+2h\alpha _0)y_0+4y_1-y_2=2h\alpha_1\\ y_{i-1}-(2+h^2q_i)y_i+y_{i+1}=h^2f_i,&i=1,2,\cdots,n-1,\\ y_{n-2}-4y_{n-1}+（3+2h\beta_0)y_n=2h\beta_1 \end{array}\right.$$
• 然后可以使用追赶法求解线性方程组
• 利用常微分方程的只是可以求出微分方程的解析解是 $$y=-\frac{1}{12}x^4-\frac{1}{6}x+1$$ 最后可以敬爱嗯差分方程的的解和真实解的图像作对比。

代码

function [] = main()
%用差分法求解二阶常微分方程
N = input('Please Inupt N:')
stap = 0;
endp = 1;
h = (endp-stap)/N;
xarr = stap:h:endp;
%% 获得差分方程的系数矩阵
A = zeros(N+1);
b = ones(N+1,1);
for t = 1:N+2
    if t==1
        A(t,t) = 1;
        b(t)   = 1;
    else
        if t<=N+1
            A(t,t-1:t+1) = [1 -2 1];
            b(t) = -h*h*xarr(t)*xarr(t);
        else
            A(t,t-1:t) = [1 (1+h)/(1-h) ];
            b(t) = h/(1-h);
        end
    end
end
%% 使用追赶法求解方程组
% disp(A);
% disp(b);
y = zeros(N+2,1);
u = zeros(N+2,1);
l = zeros(N+2,1);
x = zeros(N+2,1);
u(1) = A(1,1);
y(1) = b(1);
for t = 2:N+2
    l(t) = A(t,t-1)/u(t-1);
    u(t) = A(t,t) - l(t)*A(t-1,t);
    y(t) = b(t) - l(t)*y(t-1);
end
x(N+2)  = y(N+2)/u(N+2);
for t = N+1:-1:1;
    x(t) = (y(t)-A(t,t+1)*x(t+1))/u(t);
end
xlswrite('res.xls',x,1);
plot(xarr,x(1:N+1));
% disp(x);
end

结果

• 模拟函数图像如下
• 真实函数图像如下
  可以看出两者的差距较小，因此可以认为差分算法是有效的。

第二题

题目

用差分法求解如下 Laplace 方程 $(N=400)$

$$\left\{ \begin{array}{rcl} -\Delta u=xy,&&(x,y)\in \Omega=[0,1]\times [0,1] \\u=0,&&(x,y)\in\partial\Omega \end{array} \right.$$

分析

• 不妨记 $f(x,y) = xy, (x,y) \in \Omega = [0,1]\times [0,1]$ ,所以原方程可以写做 $-\Delta u=f(x,y)$ ，因为在边界出函数的值为零，所以是第一边值条件。
• 使用五点差分格式对方程进行差分
• 首先对函数的定义域进行剖分，根据 $N$ 的大小，可以确定网格的大小， 因此宽和高都为 $\frac{1}{N}$,以此也可以确定网点的位置。最终差分法会计算出函数在差分点的函数值，从而绘制出图形。
• 五点差分法的格式为 $$u_{ij}-\frac{1}{4}(u_{i-1,y}+u_{i,j-1}+u_{i+1,j}+u_{i,j+1})=\frac{h^2}{4}f_{ij}$$
• 利用 crs 压缩行存储以及五点差分法的格式可以获得差分方程系数矩阵的压缩存储格式以及值向量
• 利用基本的迭代法求解方程获得方程的解并绘制出图形
• 使用 freefem++ 通过变分法计算原方程的解并绘制出图形。与差分法绘制的图形进行对比。

代码

主函数 main

function [] = main()
N = input('lxease Input N:');
[value, cow, row, b]=getmatrix(N);
% for m=1:(N+1)^2
%     fprintf('\n')
%     for n = 1:(N+1)^2
%         %fprintf('%f',m)
%         temp = getvalue(m,n,value,cow,row);
%         fprintf('%f\t',temp)
%     end
% end
% b;
%% 使用基本迭代法求解方程组，（ E+ T）x = b.  then x = b-t*x
res  = solve_equation(value,cow,row,b);
%% now plot the shape
[y,x]=meshgrid(1:-1/N:0);
Z =flipud( reshape(res,N+1,N+1));
mesh(x,y,Z);
end

getmatrix 函数 主要用来获得系数矩阵的压缩行存储格式

function [value, cow, row,b] = getmatrix(N)
%% 使用addm 差值节点的地址信息。
[xsta, xend, ysta, yend] = deal(0, 1, 0, 1);%多变量同时赋值
xh = (xend-xsta)/N;
yh = (yend-ysta)/N;
addm = {};
for m = 1:N+1
    for n = 1:N+1
        xtemp = xsta + (m-1)*xh;
        ytemp = ysta + (n-1)*yh;
        addm(m,n) = {[xtemp, ytemp, xtemp*ytemp ],};
    end
end
%disp(addm);
%% 使用crs存储方法记录系数矩阵A
valuelen = (N*3-2)*N + N*(N-1)*2;
value = ones(valuelen,1);
cow = ones(N*N,1);
row = ones(valuelen,1);
tempadr = 1;%利用 tempadr 记录 value 数值变化的位置;
cowadr = 1;%利用 cowadr 记录 cow 列号变化的位置；
rowadr = 1 ;%利用 row 记录行变化的位置
b = zeros(N*N+2*N+1,1);
badr = 1;
%% 从第一行开始逐个记录系数矩阵每行的
for t = 1:N+1    
    if t==1
        for n = 1:N+1
            value(tempadr) = 1;
            tempadr = tempadr + 1;
            cow(cowadr) = tempadr;
            cowadr = cowadr +1;
            row(rowadr) = (t-1)*(N-1)+n;
            rowadr = rowadr + 1;
            b(badr) = 0;
            badr = badr + 1;
        end        
    else
        if t<=N               
            value(tempadr) = 1;
            tempadr = tempadr + 1;
            cow(cowadr) = tempadr;
            cowadr = cowadr +1;
            row(rowadr) = (t-1)*(N+1)+1;
            rowadr = rowadr + 1;
            b(badr) = 0;
            badr = badr + 1;
            for n=2:N
                value(tempadr:tempadr+4) = [-1/4  -1/4 1 -1/4 -1/4];
                tempadr = tempadr +5;
                cow(cowadr) = tempadr;
                cowadr = cowadr +1;
                row(rowadr:rowadr+4) = [(t-2)*(N+1)+n, (t-1)*(N+1)+n-1, (t-1)*(N+1)+n, (t-1)*(N+1)+n+1, t*(N+1)+n];
                rowadr = rowadr + 5;
                b(badr) = xh*xh*addm{t,n}(3);
                badr = badr + 1;
            end            
            value(tempadr) = 1;
            tempadr = tempadr + 1;
            cow(cowadr) = tempadr;
            cowadr = cowadr +1;
            row(rowadr) = t*(N+1);
            rowadr = rowadr + 1;
            b(badr) = 0;
            badr = badr + 1;                         
        else
            for n = 1:N+1
                value(tempadr) = 1;
                tempadr = tempadr + 1;
                cow(cowadr) = tempadr;
                cowadr = cowadr +1;
                row(rowadr) = (t-1)*(N+1)+n;
                rowadr = rowadr + 1;
                b(badr) = 0;
                badr = badr + 1;
            end             
        end
    end
end
end

solve_equation 主要使用基本迭代法求解方程组

function [ x ] = solve_equation( value, cow, row, b )
%This function is used to solve equation
%% get size of x;
lenth = length(b);
x0 = ones(lenth,1);
x = x0;
error = 1;
%% 方程的解为x
while error>=1e-5    
    for m=1:lenth
        temp = 0;
        if m == 1
            for t = 1:(cow(m)-1)
                if row(t) ~= m
                    n = row(t);
                    temp = temp + x0(n)*value(t);
                end
            end
        else
            for t = cow(m-1):(cow(m)-1)
                if row(t) ~= m
                    n = row(t);
                    temp = temp + x0(n)*value(t);
                end
            end
        end
        temp = b(m)- temp;
        x(m) = temp;       
    end  
   error = norm(x0-x,2);
   x0 = x;
end
end

getvalue 函数 检查系数矩阵是否正确

function[val] = getvalue(m,n,value,cow,row)
% 这个函数通过 crs 存储的 value , cow ,row 计算出矩阵第m行 n列的元素的值并返回
%%
lenth = length(cow);
if m>lenth || m<=0
    fprintf('Larger than the matrix')
else
    val = 0;
    if m==1
        for t=1:(cow(m)-1)
            if n == row(t)
                val =value(t);
                return;
            end
        end
    else
        val = 0;
        for t=cow(m-1):(cow(m)-1)
            if n == row(t)
                val =value(t);
                return;
            end
        end
    end
end

使用 freefem++ 通过变分法计算，代码如下：

real theta=4.*pi/3;
real a=2. ,b=1.;
func z=0;
border a0(t=0,1){x=1;  y=t;}
border a1(t=0,1){x=1-t;y=1;}
border a2(t=0,1){x=0;  y=1-t;}
border a3(t=0,1){x=t;y=0;}
mesh Th = buildmesh(a0(50)+a1(50)+a2(50)+a3(50));
fespace Vh(Th,P2);
Vh phi,w,f=x*y;
solve Laplace(phi,w)=int2d(Th)(dx(phi)*dx(w) + dy(phi)*dy(w))
- int2d(Th)(f*w)+on(a0,phi=z)+on(a1,phi=z)+on(a2,phi=z)+on(a3,phi=z) ;
plot(phi,wait=true);
plot(Th,wait=true, ps="membraneTh.eps");
savemesh(Th,"trace.msh");

结果

• 首先检查系数矩阵的正确性。输入的n 为2 矩阵的规模为 $9\times9$
• 程序运行结果所绘制的图形。
• 使用 freefem++ 绘制图形如下

第四次上机作业

题目

分析

代码

结果

在此输入正文