DHU Matlab Experiment【2】作业记录_第三章、第四章

写在前面

本博客用于记录(或者说是用来备份)我在2021a 高等数学实验课 梁志勇老师的课上写的程序习题
所有的代码均通过编译,matlab版本为R2016a
课程链接:高等数学实验
由于题量的关系和对于篇幅的考虑,故一章到两章会做一篇博文,更多后续章节的答案可以点击我的头像查看或者点击下面的链接查看:
第一、第二章
第五、第六章
第七章
网上参考答案

3.2.3 矩阵代数

  1. 设a=(1,2,3),b=(2,4,3), 分别计算a./b, a.\b, a/b, a\b, 分析结果的意义。
    在这里插入图片描述

  2. 求下列矩阵的行列式、逆、特征值和特征向量
    在这里插入图片描述
    (4) n阶三对角方阵
    在这里插入图片描述
    阶数n 分别为5, 50.
    注意: 由于输出结果太多, 答题框里放不下. 请大家只写出n 分别为5, 50对应的程序, 输出结果只写n=5的.
    在这里插入图片描述
    生成(对角线上元素相同的)三对角矩阵

clear all;clc;
n = 5;%当n=50时,替换此即可
A = diag(linspace(6,6,n-1),1)+...
    diag(linspace(5,5,n))+...
    diag(linspace(1,1,n-1),-1)%构造n阶三对角方阵
inv(A)%逆矩阵
[v,d] = eig(A)

在这里插入图片描述

  1. 判断上题第(1)小题是否为正定矩阵。
    判断矩阵是正定、半正定还是负定
    从参考的代码中获取我们需要的判定部分~,如下
A = [4 1 -1;3 2 -6;1 -5 3];
if issymmetric(A) % 检查矩阵是否对称
    d = eig(A) % 计算矩阵特征值
    if all(d > 0)
        disp('矩阵正定');
    end
else
    disp('矩阵不正定');
end

3.3 线性方程组的通解

  1. 用矩阵除法解下列线性方程组,并判断解的意义. 用矩阵乘法验算,只要求做第(1)、(4)小题。
    在这里插入图片描述
    (1)
    在这里插入图片描述
    注意,该题中如果用较为准确的判断,如A*x-bA*x==b将会发现在1e^-15处二者不相等,而由于’=='的判断精度是16位,因此也会报错。这个问题应该来源于精度损失,即x的有效数字为16位以上,导致matlab没有完整的地记录下该数据,导致之后的逆运算中出现误差。
    (4)
    在这里插入图片描述

  2. 求上题第(4)小题的通解
    【法一】用rref化为行最简形以后求解
    在这里插入图片描述
    进一步计算:
    在这里插入图片描述
    【法二】使用除法得到一特解,再用null得到空间的一个正交规范基
    在这里插入图片描述
    易知结果为kx+x0,与法一自洽

  3. 求下列向量组的秩和它的一个最大线性无关组,并将其余向量用该最大无关组线性表示。
    a1= (4, -3, 1,3), a2= (2, -1, 3, 5), a3= (1, -1, -1, -1), a4= (3, -2, 3, 4), a5= (7, -6, -7, 0)
    在这里插入图片描述
    将其余向量用该最大无关组线性表示:
    在这里插入图片描述

3.4 投入产出分析和基因遗传

  1. (人口流动趋势)对城乡人口流动作年度调查,发现有一个稳定的朝向城镇流动的趋势,每年农村居民的5%移居城镇而城镇居民的1%迁出,现在总人口的20%位于城镇。假如城乡总人口保持不变,并且人口流动的这种趋势继续下去,那么
    (1)一年以后住在城镇人口所占比例是多少?两年以后呢?十年以后呢?
    (2)很多年以后呢?
    (3)如果现在总人口70%位于城镇,很多年以后城镇人口所占比例是多少?
    (4)计算转移矩阵的最大特征值及对应的特征向量,与问题(2)(3)有何关系?

(1)

clear all;clc;
x0 = [0.2;0.8];%城市和农村的人口占比
a = [0.99 0.05;0.01 0.95];
x = x0;
for i = 1:1:10
    x = a*x;
    if (i== 1)|(i==2)|(i==10)
        i
        x
    end
end

高级版两个池子互相灌水的程序设计题hhh
编写的思路是找一个流动矩阵,能够迭代地进行满足题目要求的人口流动(因为矩阵的乘法规则,所以一定得这么乘)
在这里插入图片描述
检验代码结果的正确性的方法为:
因为人口没有流失,所以要将x中的元素相加,恒为1
即0.2+0.8=1;
一年后0.238+0.762=1;
两年后0.2737+0.7263=1;
十年后0.4922+0.5078=1

运行结果:
在这里插入图片描述

即一年以后城市人口占23.8%;两年后占27.37%;十年后占49.22%
(2)观察x(1)的趋势,易知其单调递增,故用增幅来考察多年后的人口占比极限

clear all;clc;
x0 = [0.2;0.8];%城市和农村的人口占比
a = [0.99 0.05;0.01 0.95];
i = 1;
x1 = a*x0;%一年后城市与农村的人口占比
while abs(x1(1)-x0(1)) > 1e-6,
    i = i+1;
    x0 = x1;%上一年的人口占比
    x1 = a*x1;%第i年后的人口占比
end
i%该年份后,增幅一直小于1e-6
x1%该年份的比例

运行结果:
在这里插入图片描述
即多年以后,城市人口占总人口的比值维持在了83.33%
(3)将x0换为题设后,换一种思路:直接查看1000年以后的城镇人口所占比例

clear all;clc;
x0 = [0.7;0.3];%城市和农村的人口占比
a = [0.99 0.05;0.01 0.95];
x = x0;
for i = 1:1:1000
    x = a*x;
end
x

运行结果:
在这里插入图片描述
即多年以后,城市人口占总人口的比值仍维持在了83.33%
将(2)中的x0修改后验证,自洽
(4)以(3)中最后得到的x为例
在这里插入图片描述
用(2)中得到的x0,进行检验,自洽

  1. (经济预测)在某经济年度内,各经济部门的投入产出表如下表3.5(单位:亿元)
    在这里插入图片描述
    假设某经济年度工业,农业及第三产业的最后需求均为17亿元,预测该经济年度工业,农业及第三产业的产出(提示:对于一个特定的经济系统而言,直接消耗矩阵和Leontief矩阵可视作不变)。
clear all;clc;%关于题设中的表格:已知B1,x1,D1求直接消耗矩阵C和Leontief矩阵A
B1 = [6 2 1;2.25 1 0.2;3 0.2 1.8];x1 = [25 5 20];
X1 = x1';D1 = [16;1.55;15];
C = B1/diag(x1);%B=C*diag(x1),右除
A = eye(3) - C;%A=E-C,E为单位矩阵
%D1 = A*X1可用于验算,但注意用D1/X1不能求得A
%因为X1不为方阵,不能求其逆(见P49式子3.4上方对于可逆的充要条件)
%下求对某经济年度的预测
D = [17;17;17];
X = A\D%已用rank(A)=rank([A,D])=n证明X为唯一解

运行结果:
在这里插入图片描述
即该经济年度工业的总产值为37.5696;农业的总产值为25.7862;第三产业的总产值为24.7690

4.2 函数零点、极值和最小二乘拟合

  1. 求下列多项式的所有根,并进行验算。(请按题号标注习题解答,只用做第(2)、(3)小题)
    在这里插入图片描述
    (2)
    在这里插入图片描述
    (3)
    在这里插入图片描述
    在这里插入图片描述
  2. 求方程
    在这里插入图片描述
    的正根
clear all;clc;
fun = @(x)x.*log(sqrt(x.^2-1)+x)-sqrt(x.^2-1)-0.5*x;
fplot(fun,[1,3]);%定义域x^2-1>0
grid on;
%可见上述方程有唯一正根在2.1附近,取x0=2.1
[x,f,h] = fsolve(fun,2.1)

运行结果:

在这里插入图片描述

  1. 求解下列非线性方程组在原点附近的根
    在这里插入图片描述
    函数:文件名命名为“ex3_fun.m”
function f = fun(x)
%x(1):x    x(2):y    x(3):z
f(1) = 9*x(1).^2+36*x(2)^2+4*x(3).^2-36;
f(2) = x(1).^2-2*x(2)^2-20*x(3);
f(3) = 16*x(1)-x(1).^3-2*x(2)^2-16*x(3).^2;
end

主程序

clear;[x,f,h]=fsolve(@ex3_fun,[0,0,0])

运行结果:
在这里插入图片描述
易知得到的解为x=0.1342,y=0.9972,z=-0.0985,根据误差且h>0知结果可靠

  1. 作出下列函数图形,观察所有的局部极大, 局部极小和全局最大, 全局最小值点的粗略位置; 并用MATLAB函数fminbnd和fminsearch求各极值点的确切位置
    在这里插入图片描述
clear;clc;
fun = @(x)x.^2.*sin(x.^2-x-2);
fplot(fun,[-2,2]);grid on;
%观察得局部极大的粗略位置x0=-1.5,0
%局部极小粗略位置x0=-2,-0.5,1.5
%全局最大粗略位置x0=-1.5
%全局最小的粗略位置x0=-2
[x(1),f(1)] = fminbnd(fun,-2,-1.5);%局部极小值
[x(2),f(2)] = fminsearch(fun,-0.7);%局部极小值
[x(3),f(3)] = fminsearch(fun,1.5);%局部极小值
fun2 = @(x)-x.^2.*sin(x.^2-x-2);
[x(4),f(4)] = fminbnd(fun,-1.7,-1.2);%局部极大值
[x(5),f(5)] = fminsearch(fun,0);%局部极大值
%输出部分
disp([x(1),x(2),x(3)]);%局部极小值
disp([x(4),x(5)]);%局部极小值
x(find(f==min(f)))%全局最小值
x(find(f==max(f)))%全局最大值

运行结果:
在这里插入图片描述
在这里插入图片描述

  1. 考虑函数
    在这里插入图片描述
    (1)作出f(x,y)在-2<x<1, -7<y<1 的图,观察极值点的位置;
    (2) 用MATLAB函数fminsearch求极值点和极值。
clear all;clc;
x = -2:0.1:1;y = -7:0.1:1;
[x,y] = meshgrid(x,y);
z = y.^3/9+3.*x.^2.*y+9.*x.^2+x.*y+9+y.^2;
mesh(x,y,z);grid on;
xlabel('x轴');ylabel('y轴');zlabel('z轴');
fun = @(x) x(2).^3/9+3.*x(1).^2.*x(2)+9.*x(1).^2+x(1).*x(2)+9+x(2).^2;
[x1 f1] = fminsearch(fun,[-1.9,-6.9]);%极小值
fun1 = @(x) -(x(2).^3/9+3.*x(1).^2.*x(2)+9.*x(1).^2+x(1).*x(2)+9+x(2).^2);
[x2 f2] = fminsearch(fun,[0,-5]);%极大值
x1%极小值
x2%极大值

就算做了图,因为是三维的,也看不太清极值点
在这里插入图片描述
在这里插入图片描述

4.3.3 迭代法

  1. 假定某天的气温变化记录如第二章习题5,试用最小二乘方法找出这一天的气温变化规律。考虑下列类型函数, 作图比较效果,并比较误差平方和的大小。
    在这里插入图片描述
    (1) 二次函数;
    (2) 三次函数;
    (3) 钟形函数在这里插入图片描述

(1)

clear;clc;
x = 0:1:24;y = [15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];
plot(x,y,'ro');xlabel('时刻t(h)');ylabel('温度(℃)');hold on;
c1 = polyfit(x,y,2);%二次函数多项式拟合
fun_out1 = @(x) c1(1)*x.^2+c1(2)*x+c1(3)
fplot(fun_out1,[0,24],'b');%作图比较
sum((feval(fun_out1,x)-y).^2)%计算误差平方和

在这里插入图片描述
在这里插入图片描述
(2)

clear;clc;
x = 0:1:24;y = [15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];
plot(x,y,'ro');xlabel('时刻t(h)');ylabel('温度(℃)');hold on;
c2 = polyfit(x,y,3);%二次函数多项式拟合
fun_out2 = @(x) c2(1)*x.^3+c2(2)*x.^2+c2(3)*x+c2(4)
fplot(fun_out2,[0,24],'g');%作图比较
sum((feval(fun_out2,x)-y).^2)%计算误差平方和

在这里插入图片描述
在这里插入图片描述
(3)

clear;clc;
x = 0:1:24;y = [15 14 14 14 14 15 16 18 20 22 23 25 28 31 32 31 29 27 25 24 22 20 18 17 16];
plot(x,y,'ro');xlabel('时刻t(h)');ylabel('温度(℃)');hold on;
fun = @(c,x) c(1).*exp(c(2)*(x-14).^2);
[c3,Q] = lsqcurvefit(fun,[0,0],x,y)
fun_out3 =  @(x) c3(1).*exp(c3(2)*(x-14).^2);
fplot(fun_out3,[0,24]);hold off;

在这里插入图片描述
在这里插入图片描述
比较三种拟合方式:
在这里插入图片描述
其中,红色圆圈代表原数据;蓝色线代表(1);绿色代表(2);黑色代表(3);
从图像上看,使用法(3),即用钟形函数,拟合得更接近原曲线;从误差平方和的角度分析,法(2),即用三次函数拟合得更接近原曲线

  1. (弦截法)牛顿迭代法是一种速度很快的迭代方法,但是它需要预先求得导函数。若用差商代替导数,可得下列弦截法
    在这里插入图片描述
    这一迭代法需要两个初值x0, x1,编写一个通用的弦截法计算机程序并用以解习题2。
    习题2:
    在这里插入图片描述
clear;clc;
x0 = 1;%函数的定义域的起始位置
x1 = x0 + 0.01;
fun = inline('x.*log(sqrt(x.^2-1)+x)-sqrt(x.^2-1)-0.5.*x','x');
fx0 = feval(fun,x0);
fx1 = feval(fun,x1);
while abs(fx0) > 1e-6,
    x2 = x1 - (x1-x0)/(fx1-fx0)*fx1;%xk+1
    %更新参数以进入下一次循环
    x0 = x1;%xk
    x1 = x2;%xk+1
    fx0 = feval(fun,x0);%f(xk)
    fx1 = feval(fun,x1);%f(xk+1)
end
x0

运行结果:
在这里插入图片描述

与习题2中的结果比较:
在这里插入图片描述
答案一致

4.4.5 购房贷款的利率

  1. (化学反应平衡) 一等克分子数一氧化碳(CO)和氧气(O2)的混合物在300K和5bar压力下达到平衡,理论反应方程式为
    在这里插入图片描述
clear;clc;
fun = @(x)(1-x).*sqrt(10.52+x)./(x.*sqrt(1+x).*5-3.06);
fplot(fun,[0,1]);hold off;%观察解的大致位置
a1 = fzero(fun,0.5)
feval(fun,a1)%验算

运行结果:
在这里插入图片描述

  1. (月还款额)作为房产公司的代理人,你要迅速准确回答客户各方面的问题。现在有个客户看中了你公司一套建筑面积为180平方米,每平方单价7500元的房子。他计划首付30%,其余70%用20年按揭贷款(贷款年利率5.04%)。请你提供下列信息:房屋总价格、首付款额、月付还款额。如果其中10万元为公积金贷款(贷款年利率4.05%)呢?
clear;clc;
x = 180 * 7500/10000%房屋总价格,单位:万元
x0_1 = x * 0.3%首付款额
x0 = x * 0.7;%总还款额
r = 0.0504/12;%月利率
N = 20*12;%贷款期限(月)
a = (1+r)^N*r*x0/((1+r)^N-1)%单位:万元
%若其中10万元为公积金贷款(贷款年利率4.05%)
x1 = 10;%公积金贷款
x2 = x0-x1;%总还款额
r1 = 4.05/100/12;
a1 = (1+r1)^N*r1*x1/((1+r1)^N-1);
a2 = (1+r)^N*r*x2/((1+r)^N-1);
a0 =a1+a2%用一部分公积金贷款的月还贷

运行结果:
在这里插入图片描述
即房屋总价格:135万元;首付款额:40.5万元;月付还款额:6257元
若其中10万元为公积金贷款(贷款年利率4.05%),月付还款额:6204元

  1. (栓牛鼻的绳子)农夫老李有一个半径10米的圆形牛栏,里面长满了草,老李要将家里一头牛栓在一根栏桩上,但只让牛吃到一半草,他想让上大学的儿子告诉他,栓牛鼻的绳子应为多长?
    思路:
    两圆的示意图如下:
    在这里插入图片描述
    求蓝色部分的面积,用两个扇形-四边形O1AO2B,即两个扇形-两个三角形O1AO2即可
    (按题设画的图比较难看出来面积关系)
    在这里插入图片描述
    费了老大劲用角度制算的,忘记matlab里三角函数用的是弧度制了,后悔
    但至少说明了用角度制也能算…稍微改一改也能送进matlab里运行,就是一开始忘记三角函数用弧度制的时候debug用了不少时间
clear;clc;
r = 10;%单位:米
S = @(th)-pi*(180-th)/180.*cos(th*pi/180)-sin(th*pi/180)+1/2*pi;
th = fsolve(S,60)%theta粗略范围:0-90
R = sqrt(2*r^2-2*r^2*cos(th*pi/180))%R的粗略范围10-20

注意设不同的角获得的式子不同,但是做出来的R都是一致的
运行结果
在这里插入图片描述

  • 52
    点赞
  • 163
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 2
    评论

“相关推荐”对你有帮助么?

  • 非常没帮助
  • 没帮助
  • 一般
  • 有帮助
  • 非常有帮助
提交
评论 2
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

鱼犬

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值