中国剩余定理matlab代码
function CHN_Rem_Thm
a=input('输入数组a:');
m=input('输入数组m:');
%a=[2 8 5];
%m=[17 11 19];
if length(a)~=length(m)
error('Length of a and length of m must
agree!');
return;
end
if flg==0
error('不互素');
return;
end
M=prod(m)
u=ones(size(m));
for k=1:length(m)
u(k)=inver(M/m(k),m(k));
end
x=a*u';
x=mod(x,M)
function d=inver(m,n)
for k=0:(n-1)
if mod(k*m,n)==1
d=k*m;break;
end
end
function flg=coprime_chk(m)
flg=1;
n=length(m);
for k=1:n-1
for p=(k+1):n
if gcd(m(k),m(p))~=1
flg=0;
break;
end
end
end
数学建模常用MATLAB程序代码
V1到V11的最短路径及其长度：P53>> a(2,3)=6;a(2,5)=1;a(3,4)=7;a(3,5)=5;a(3,6)=1;a(3,7)=2;a(4,7)=9;a(5,6)=3;a(5,8)=2;a(5,9)=9;a(6,7)=4;a(6,9)=6;a(7,9)=3;a(7,10)=1;a(8,9)=7;a(8,11)=9;a(9,10)=1;a(9,11)=2;a(10,11)=4;a=a';[i,j,v]=find(a);b=sparse(i,j,v,11,11);[x,y,z]=graphshortestpath(b,1,11,'Directed',false)x =13y =1 2 5 6 3 7 10 9 11z =0 1 6 1 2 5 3 5 10 7 9灰色预测GM(1,1)模型程序代码：P375x0=[71.1 72.4 72.4 72.1 71.4 72 71.6]';n=length(x0);lamda=x0(1:n-1)./x0(2:n);range=minmax(lamda');x1=cumsum(x0)B=[-0.5*(x1(1:n-1)+x1(2:n)),ones(n-1,1)]Y=x0(2:n);u=B\Yx=dsolve('Dx+a*x=b','x(0)=x0')x=subs(x,{'a','b','x0'},{u(1),u(2),x0(1)})yuce1=subs(x,'t',[0:n-1])y=vpa(x,n-1)yuce=[x0(1),diff(yuce1)]epsilon=x0'-yucedelta=abs(epsilon./x0')rho=1-(1-0.5*u(1))/(1+0.5*u(1))*lamda'x1 =71.1143.5215.9288359.4431.4503B =-107.3 1-179.7 1-251.95 1-323.7 1-395.4 1-467.2 1u =0.002343872.657x =(b - exp(-a*t)*(b - a*x0))/ax =41884064296049311744/1351100916648171 - (41788001020875628544*exp(-(1351100916648171*t)/576460752303423488))/135110091664 8171yuce1 =71.1 143.51 215.74 287.81 359.71 431.44 503y =31000.0 - 30928.9*exp(-0.00234379*t)yuce =71.1 72.406 72.236 72.067 71.898 71.73 71.562epsilon =0 -0.0057414 0.16376 0.032871 -0.49842 0.2699 0.037824delta =0 7.9302e-05 0.0022619 0.00045592 0.0069806 0.0037486 0.00052827rho =0.020255 0.002341 -0.0018101 -0.0074399 0.010655 -0.0032325用Krusksl算法求最小生成树的Matlab程序代码：P47-48a(1,2)=50;a(1,3)=60;a(2,4)=65;a(2,5)=40;a(3,4)=52;a(3,7)=45;a(4,5)=50;a(4,6)=30;a(4,7)=42;a(5,6)=70;>> [i,j,b]=find(a);>> data=[i';j';b'];index=data(1:2,:);>> loop=length(a)-1;>> result=[];>> while length(result)<looptemp=min(data(3,:));flag=find(data(3,:)==temp);flag=flag(1);v1=index(1,flag);v2=index(2,flag);if v1~=v2result=[result,data(:,flag)];endindex(find(index==v2))=v1;data(:,flag)=[];index(:,flag)=[];end>> resultresult =6 57 7 2 54 2 4 3 1 430 40 42 45 50 50求得最小生成树的边集为{v6v4, v5v2, v7v4, v7v3, v2v1, v5v4}。
matlab数学建模程序代码
matlab数学建模程序代码
当进行数学建模时，MATLAB是一个强大的工具，用于实现和测试模型。
下面是一个简单的MATLAB代码示例，演示如何使用MATLAB进行一维线性回归建模：
```matlab
%生成示例数据
x=[1,2,3,4,5];
y=[2.8,3.9,4.8,5.5,6.3];
%进行一维线性回归
coefficients=polyfit(x,y,1);
slope=coefficients(1);
intercept=coefficients(2);
%绘制原始数据和回归线
scatter(x,y,'o','DisplayName','原始数据');
hold on;
plot(x,polyval(coefficients,x),'r-','DisplayName','回归线');
hold off;
%添加标签和图例
xlabel('X轴');
ylabel('Y轴');
title('一维线性回归建模示例');
legend('show');
%输出回归方程的系数
fprintf('回归方程：y=%.2fx+%.2f\n',slope,intercept);
```
此代码生成了一些示例数据，然后使用一维线性回归对数据进行建模。
回归方程的系数将被计算，并且原始数据与回归线将在图上显示。
请注意，这只是一个简单的示例，实际上，你可能需要根据你的具体问题修改代码。
中国余数定理c语言
中国余数定理c语言摘要：1.引言2.中国余数定理简介3.中国余数定理在C 语言中的应用4.总结正文：1.引言中国余数定理，又称孙子定理，是我国古代数学家孙子所提出的一个关于同余方程组的定理。
随着现代计算机科学的飞速发展，中国余数定理在计算机编程领域也有着广泛的应用，特别是在C 语言中。
本文将简要介绍中国余数定理，以及如何在C 语言中运用这一定理。
2.中国余数定理简介中国余数定理是一种解决同余方程组的方法，它可以求解一组同余方程在模某个数时的非负整数解。
该定理的表述如下：设a_1, a_2, ..., a_n 两两互质，m 是a_1, a_2, ..., a_n 的最大公约数，那么同余方程组：x ≡ a_1 (mod m)x ≡ a_2 (mod m)...x ≡ a_n (mod m)有唯一解x ≡ 0 (mod m)。
3.中国余数定理在C 语言中的应用在C 语言中，中国余数定理可以应用于多种计算场景，例如大数运算、模运算等。
下面以一个简单的示例来说明如何在C 语言中使用中国余数定理求解同余方程组：```c#include <stdio.h>#include <stdbool.h>bool chinese_remainder(int a[], int m[], int n, int M) {int x = 0;for (int i = 0; i < n; ++i) {if (m[i] == 1) {x = (x + a[i]) % M;} else {int u = chinese_remainder(a, m + 1, n - 1, M / m[i]);x = (x * m[i] + u) % M;}}return x == 0;}int main() {int a[] = {2, 3, 5};int m[] = {5, 7, 11};int n = sizeof(a) / sizeof(a[0]);int M = 1;for (int i = 0; i < n; ++i) {M *= m[i];}if (chinese_remainder(a, m, n, M)) {printf("同余方程组有解");} else {printf("同余方程组无解");}return 0;}```在这个示例中，我们定义了一个名为`chinese_remainder` 的函数，用于求解给定的同余方程组。
matlab代码
n=7:12;y=1./abs(n-6);plo t(n,y,'r*','Mar kersi ze',20)gr id on//画图x=l inspa ce(0,2*pi,100);y=si n(x);plot(x,y);加注解x=l inspa ce(0,2*pi,100);>> y=sin(x);>> plo t(x,y);>> x=li nspac e(0,2*pi,100);>> y=sin(x);>> plot(x,y,'r',x,y,'*');>> leg end('y=sin(x)')>> t ext(1.4,0.91,'波峰')x label('Inp ut Va lue');>>ylabe l('Fu nctio n Val ue');>> t itle('一个正弦函数');//a xis命令控制坐标轴x=l inspa ce(0,2*pi,100);y=si n(x);z=co s(x);plot(x,y,x,z);axis on>> axi s('sq uare')>>axis('equa l')//plo t函数t=(0:p i/100:pi);y1=s in(t)*[1,-1];y2=sin(t).*sin(9*t);t3=pi*(0:9)/9;y3=si n(t3).*sin(9*t3);pl ot(t,y1,'r:',t,y2,'b',t3,y3,'b o')a xis([0,pi,-1,1])//hold on 命令不会取消之前画的曲线>> hold on>> y4=-1:0.2:1;>> pl ot(1.5,y4,'*')//同一窗口用不同的坐标系绘制几组数据subp lot(2,2,1);plot(x,si n(x));>>subpl ot(2,2,2);plot(x,cos(x));>> s ubplo t(2,2,3);p lot(x,sinh(x));>> s ubplo t(2,2,4);p lot(x,cosh(x));//f igure函数x=lins pace(0,2*p i,30);>>figur e>>plot(x,sin(x));>> x label('Inp ut Va lue');>>ylabe l('Fu nctio n Val ue');>> t itle('一个正弦函数');>> l egend('y=s in(x)')>> figu re>> plot(x,si nh(x));>> clos e//直方图close all>> x=1:10;>> y=rand(size(x));>> b ar(x,y);>> x=1:10x=lin space(0,2,20)*p iy=s in(x)e=st d(y)*ones(size(x));error bar(x,y,e)//柄图x=l inspa ce(0,10,50);>> y=si n(x).*exp(-x/3);>>stem(x,y)//阶梯图> x=lin space(0,10,50);y=si n(x).*exp(-x/3);>>stair s(x,y);//饼图x=[11.4,23.5,35.4,15.6];>> ex plode=zero s(siz e(x));>>[c,of fset]=min(x);>> exp lode(offse t)=1;>> p ie(x,explo de);//频数累计柱状图x=ra ndn(5000,1);>> hist(x,20);//矢量图//羽状图thet a=90:-10:0;>>r=one s(siz e(the ta));>> [u,v]=pol2c art(t heta*pi/180,r*10);>> fea ther(u,v)>> ax is eq ual//极坐标图th eta=l inspa ce(0,2*pi);>>r=cos(4*th eta);>> p olar(theta,r);//精确绘图fp lot('sin(1./x)',[0.02,0.2]); 0.02,0.2是绘图范围三维绘图三维线性图形>> t=(0:0.02:2)*pi;>> x=sin(t);>> y=c os(t);>>z=cos(2*t);>>plot3(x,y,z,'b-',x,y,z,'b d');>> vi ew([-82,58]);>> box on;>> le gend('链','宝石');//三维曲面x=-4:0.2:4;>>y=x;>> [X,Y]=m eshgr id(x,y);>> Z=X.^2/9+Y.^2/9;>> sub plot(1,2,1);>> mesh(X,Y,Z)>> titl e('椭圆抛物面网线图')>> sub plot(1,2,2);>> surf(X,Y,Z)>> tit le('椭圆抛物面网面图')//设置视角vi ew(az,el):az:方位角，el:仰视角v iew([x,y,z]):由坐标指定观察点，即确定视角。
中国剩余定理(孙子问题)PPT课件( 13页)
•
8、不要活在别人眼中，更不要活在别人嘴中。世界不会因为你的抱怨不满而为你改变，你能做到的只有改变你自己！
•
9、欲戴王冠，必承其重。哪有什么好命天赐，不都是一路披荆斩棘才换来的。
•
10、放手如拔牙。牙被拔掉的那一刻，你会觉得解脱。但舌头总会不由自主地往那个空空的牙洞里舔，一天数次。不痛了不代表你能完全无视，留下的那个空缺永远都在，偶尔甚至会异常挂念。适应是需要时间的，但牙总是要拔，因为太痛，所以终归还是要放手，随它去。
•
6、人性本善，纯如清溪流水凝露莹烁。欲望与情绪如风沙袭扰，把原本如天空旷蔚蓝的心蒙蔽。但我知道，每个人的心灵深处，不管乌云密布还是阴淤苍茫，但依然有一道彩虹，亮丽于心中某处。
•
7、每个人的心里，都藏着一个了不起的自己，只要你不颓废，不消极，一直悄悄酝酿着乐观，培养着豁达，坚持着善良，只要在路上，就没有到达不了的远方！
•
14、给自己一份坚强，擦干眼泪；给自己一份自信，不卑不亢；给自己一份洒脱，悠然前行。轻轻品，静静藏。为了看阳光，我来到这世上；为了与阳光同行，我笑对忧伤。
•
15、总不能流血就喊痛，怕黑就开灯，想念就联系，疲惫就放空，被孤立就讨好，脆弱就想家，不要被现在而蒙蔽双眼，终究是要长大，最漆黑的那段路终要自己走完。
引入记号：m被3除余2用符号表示为Mod(m，3)
＝2；m被5除余3用符号表示为Mod(m，5)＝3；m被 7除余3用符号表示为Mod(m，7)＝2
流程图
伪代码
m2 While Mod (m,3)≠2
or Mod (m,5)≠3 or Mod (m,7)≠2 m m＋1 End While Print m
•
14、一个人的知识，通过学习可以得到;一个人的成长，就必须通过磨练。若是自己没有尽力，就没有资格批评别人不用心。开口抱怨很容易，但是闭嘴努力的人更加值得尊敬。
MATLAB软件应用 解题代码 线性代数
2%求全部4位水仙花数m=1000:9999%生成所有四位数的m向量m1=rem(m,10)%求m除以10的余数赋给m1m2=rem(fix(m/10),10)%先调用fix，对m/10的结果取整，再求取整后的数/10的余数m3=rem(fix(m/100),10)%先调用fix，对m/100的结果取整，再求取整后的数/10的余数m4=fix(m/1000)%对m/1000取整赋予m4k=find(m==m1.*m1.*m1.*m1+m2.*m2.*m2.*m2+m3.*m3.*m3.*m3+m4.*m4.*m4.*m4) %调用find函数，在向量m中找到水仙花数的序号赋给变量ks=m(k)%输出水仙花数%反序x='Kat is beautiful'%生成字符串result=strrep(x,'is',x(6:-1:5))%调用strrep函数将x中的子字符串is替换为利用冒号表达式表示的反序，初项为6，末项为5，步长为-13A=magic(7)%产生一个7阶魔方矩阵A([2 6],:)=A([6 2],:)%交换矩阵A的第一行和第二行B=A%得到矩阵Bsum(B(1,:))%求B矩阵第一行元素之和sum(B(1,:))==sum(B(2,:))%验证B矩阵各行元素之和是否相等sum(B(2,:))==sum(B(3,:))sum(B(3,:))==sum(B(4,:))sum(B(4,:))==sum(B(5,:))sum(B(5,:))==sum(B(6,:))sum(B(6,:))==sum(B(7,:))sum(B(:,1))%求B矩阵第一列元素之和sum(B(:,1))==sum(B(:,2))%验证B矩阵各列元素之和是否相等sum(B(:,2))==sum(B(:,3))sum(B(:,3))==sum(B(:,4))sum(B(:,4))==sum(B(:,5))sum(B(:,5))==sum(B(:,6))sum(B(:,6))==sum(B(:,7))sum(diag(B))%求B矩阵主对角线元素之和B1=flipud(B)%将B矩阵上下翻转，B矩阵的副对角线就变成了B1矩阵的主对角线sum(diag(B1))%求B1矩阵主对角线元素之和即求B矩阵副对角线元素之和sum(B(1,:))==sum(B(:,1))%验证B矩阵每行、每列元素之和是否相等sum(diag(B))==sum(diag(B1))%验证B矩阵主、副对角线元素之和是否相等sum(B(1,:))==sum(diag(B))%验证B矩阵每行、每列元素之以及主、副对角线元素之和是否相等B(:,[2 6])=A(:,[6 2])%交换矩阵B的第一列和第二列C=B%得到矩阵Csum(C(1,:))%求C矩阵第一行元素之和sum(C(1,:))==sum(C(2,:))%验证C矩阵各行元素之和是否相等sum(C(2,:))==sum(C(3,:))sum(C(3,:))==sum(C(4,:))sum(C(4,:))==sum(C(5,:))sum(C(5,:))==sum(C(6,:))sum(C(6,:))==sum(C(7,:))sum(C(:,1))%求C矩阵第一列元素之和sum(C(:,1))==sum(C(:,2))%验证C矩阵各列元素之和是否相等sum(C(:,2))==sum(C(:,3))sum(C(:,3))==sum(C(:,4))sum(C(:,4))==sum(C(:,5))sum(C(:,5))==sum(C(:,6))sum(C(:,6))==sum(C(:,7))sum(diag(C))%求C矩阵主对角线元素之和C1=flipud(C)%将C矩阵上下翻转，C矩阵的副对角线就变成了C1矩阵的主对角线sum(diag(C1))%求C1矩阵主对角线元素之和即求C矩阵副对角线元素之和sum(C(1,:))==sum(C(:,1))%验证C矩阵每行、每列元素之和是否相等sum(diag(C))==sum(diag(C1))%验证C矩阵主、副对角线元素之和是否相等sum(C(1,:))==sum(diag(C))%验证C矩阵每行、每列元素之以及主、副对角线元素之和是否相等4x=[2,4,8;0,-6,-4;8,1,7]%自定义一个非奇异的3阶方阵xn=norm(x) %直接求x的2-范数max(eig(x.'*x))^(1/2)%求矩阵x的转置乘矩阵x本身所得的矩阵的最大特征值的平方根%两者的结果相等x=[0,6,6,1,1,6,6,1,1,6,6,0;0,0,1,1,5,5,6,6,10,10,11,11]%生成结点坐标矩阵x，第一行为横坐标，第二行为纵坐标A1=[2,1;1,1]%定义变换矩阵A1A2=[0,-1;1,0]%定义变换矩阵A2A3=[2,0;0,2]%定义变换矩阵A3y1=A1*x%再利用A1对x进行变换，得到y1矩阵y2=A2*x%再利用A2对x进行变换，得到y2矩阵y3=A3*x%再利用A3对x进行变换，得到y3矩阵%分别绘制变换前后的图形subplot(2,2,1)fill(x(1,:),x(2,:),'r')subplot(2,2,2)fill(y1(1,:),y1(2,:),'r')subplot(2,2,3)fill(y2(1,:),y2(2,:),'r')subplot(2,2,4)fill(y3(1,:),y3(2,:),'r')5%判断n是否为质数n1 = input('请输入一个不大于100的整数：','s');n = str2num(n1);%因为默认输入的是ASCII码值，所以要有这一步，转换成数字if rem(n,2)==0 & n/2~=1 %表示如果输入的数字被2整除，则输出'您所输入的数为合数'，且把2排除disp('您所输入的数为合数')elseif rem(n,3)==0 & n/3~=1%条件表示如果输入的数字是奇数，且被3整除，则输出'您所输入的数为合数'，且把3排除disp('您所输入的数为合数')elseif rem(n,3)~=0 & rem(n,5)==0 & n/5~=1%条件表示如果输入的数字不被2、3整除，但被5这整除，则输出'您所输入的数为合数'，且把5排除disp('您所输入的数为合数')elseif rem(n,3)~=0 & rem(n,5)~=0 & rem(n,7)==0 & n/7~=1%条件表示如果输入的数字不被2、3、5整除，但被7这整除，则输出'您所输入的数为合数'，且把7排除disp('您所输入的数为合数')else %以上条件都不满足，则输出num本身disp(n)end%判断n是否为两个质数的乘积n1 = input('请输入一个不大于100的整数：','s');n = str2num(n1);%因为默认输入的是ASCII码值，所以要有这一步，转换成数字if rem(n,2)==0 & isprime(n/2) %表示如果输入的数字被2整除且除以2后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是2和',num2str(n/2),'的乘积']) %输出‘n是2和n/2的乘积’elseif rem(n,3)==0 & isprime(n/3) %表示如果输入的数字被3整除且除以3后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是3和',num2str(n/3),'的乘积']) %输出‘n是3和n/3的乘积’elseif rem(n,3)~=0 & rem(n,5)==0 & isprime(n/5) %表示如果输入的数字被5整除且除以5后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是5和',num2str(n/5),'的乘积']) %输出‘n是5和n/5的乘积’elseif rem(n,3)~=0 & rem(n,5)~=0 & rem(n,7)==0 & isprime(n/7) %表示如果输入的数字被7整除且除以7后为质数，则判定输入的数为两个质数的乘积disp([num2str(n),'是7和',num2str(n/7),'的乘积']) %输出‘n是7和n/7的乘积’else %以上条件都不满足，则输出'您输入的数不能分解为两个质数的乘积' disp('您输入的数不能分解为两个质数的乘积')end6%用辛普森法计算四分之一圆的面积，从而计算圆周率的近似值disp('以下是辛普森法求圆周率')a=0; %给a赋值为0b=1; %给b赋值为1n=input('n=?'); %输入n的值，即把区间分成n等份h=(b-a)/n; %计算每个区间的宽度，赋给变量hx=a:h:b %产生自变量向量x，x中包含了要取的n+1个自变量的值f=sqrt(1-x.*x);%得到一个与x对应的函数值向量赋给变量f,(这里一定要用点乘)s=[] %定义一个空矩阵sfor k =1:2:n %循环n次，步长为2s1=(f(k)+4*f(k+1)+f(k+2))*h/3; %求以抛物线弧段为曲边，以[Xk,Xk+2]为底的曲边梯形面积s=[s,s1]; %将s1添加到矩阵s中，循环n次后，s里面的元素分别时endpai=4*sum(s) %调用sum函数求s各元素之和，即整个曲边梯形的面积，再乘以4即为pi8figure(1)u1=linspace(-pi,pi,100);u2=linspace(-2*pi,2*pi,200);u3=linspace(0,2*pi,400);plot(u1,sin(2*u1)+1,'b-',u2,cos(u2),'r--',u3,sin(u3)+cos(2*u3),'k-.')figure(1)x=linspace(0,0.3);y=x.*sin(1./x);plot(x,y,'b')title('figure(1)')figure(2)r=1;t=0:6*pi;x1=r*(t-sin(t));y1=r*(1-cos(t));plot(x1,y1,'r')title('figure(2)')figure(3)syms x yezplot('4*x^2+9*y^2=36')grid ontitle('figure(3)')10figuresubplot(1,2,1);t=-1:0.1:1;[X,Y]=meshgrid(t);Z=X.^3+Y.^2;mesh(X,Y,Z);hold onsurf(X,zeros(size(Y)),Y.^3)title('{z=x^{3}+y^{2}}与x0z平面的相交图')subplot(1,2,2);t=-1:0.1:1;[X,Y]=meshgrid(t);Z=X.^2+Y.^5;mesh(X,Y,Z);hold onsurf(zeros(size(Y)),Y,X.^5)title('{z=x^{2}+y^{5}}与y0z平面的相交图')11%画月亮t=linspace(0,2*pi,100);subplot(2,2,1)x=sin(t);y=cos(t);p1=y>0.5;y(p1)=NaN;plot(x,y)hold onx11=sqrt(3)*sin(t)./2;y11=sqrt(3)*(0.5+cos(t))./2;p11=y11>0.5;y11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid onsubplot(2,2,2)x=sin(t);y=cos(t);p2=x>0.5;x(p2)=NaN;plot(x,y)hold onx11=sqrt(3)*(0.5+sin(t))./2;y11=sqrt(3)*cos(t)./2;p11=x11>0.5;x11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid onsubplot(2,2,3)x=sin(t);y=cos(t);p3=y<-0.5;y(p3)=NaN;plot(x,y)hold onx11=-sqrt(3)*sin(t)./2;y11=-sqrt(3)*(0.5+cos(t))./2; p11=y11<-0.5;y11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid onsubplot(2,2,4)x=sin(t);y=cos(t);p4=x<-0.5;x(p4)=NaN;plot(x,y)hold onx11=-sqrt(3)*(0.5+sin(t))./2; y11=-sqrt(3)*cos(t)./2;p11=x11<-0.5;x11(p11)=NaN;plot(x11,y11)axis([-1.1,1.1,-1.1,1.1])axis squaregrid on13% (a)x1=[1,3,5];p=poly(x1);p=p*3;px=poly2str(p,'x')x2=[2,4];q=poly(x2);q=q*-2;qx=poly2str(q,'x')% (b)A=[2,3;1,1];B=[2,1;1,1];pA=polyvalm(p,A)qA=polyvalm(q,A)pB=polyvalm(p,B)qB=polyvalm(q,B)%点乘的话，p(A)q(A)=q(B)p(A) pA.*qBqB.*pA%乘的话，p(A)q(A)和q(B)p(A)不相等pA*qBqB*pA% (c)a=0;b=10;f=diff(polyval(polyint(p),[a b]))p=[3,-27,69,-45]px=poly2str(p,'x')q=[-2,12,-16]qx=poly2str(q,'x')r=[1,1]rx=poly2str(r,'x')q1=[0,0,q]f=conv(p,r)+q1x=roots(f)14%算积分%【a】fun1=@(x) x+x.^2+x.^3;I1=rectangular(fun1,-1,1,1000)x=linspace(-1, 1, 100);y=x+x.^2+x.^3;I2=trapz(x,y)I3=quad(fun1,-1,1)%【b】fun2=@(x,y) sin(y.*(x+y)./(x.^2+4)) I=quad2d(fun2,1,10,1,10)function f=rectangular(fun,a,b,n)%矩形法:h=(b-a)/n;x=a:h:b;y=x;for i=2:n+1y(i)=fun((x(i)+x(i-1))/2);endf=h*sum(y(1:end));endsyms x y zeq1 = 3*x + 2*y - z == 4;eq2 = x - y + z == 1;eq3 = -x + 4*y + 5*z == 8;[x,y,z] = solve([eq1,eq2,eq3], [x,y,z])14%算积分%【a】fun1=@(x) x+x.^2+x.^3;I1=rectangular(fun1,-1,1,1000)x=linspace(-1, 1, 100);y=x+x.^2+x.^3;I2=trapz(x,y)I3=quad(fun1,-1,1)%【b】fun2=@(x,y) sin(y.*(x+y)./(x.^2+4)) I=quad2d(fun2,1,10,1,10)function f=rectangular(fun,a,b,n)%矩形法:h=(b-a)/n;x=a:h:b;y=x;for i=2:n+1y(i)=fun((x(i)+x(i-1))/2);endf=h*sum(y(1:end));endsyms x y zeq1 = 3*x + 2*y - z == 4;eq2 = x - y + z == 1;eq3 = -x + 4*y + 5*z == 8;[x,y,z] = solve([eq1,eq2,eq3], [x,y,z])。
matlab程序代码及学习小结
1.信号傅里叶分解代码：（1）离散傅里叶变换N=256;dt=0.02;%数据的个数和采样间隔n=0:N-1; t=n*dt;%序号序列和时间序列x=sin(2*pi*t);%合成信号m=floor(N/2)+1;%分解a，b的最大序号值为分解的N/2个参数加参数a0a=zeros(1,m);b=zeros(1,m);%产生a，b两个为零的序列for k=0:m-1,for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N);b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);%MATLAB中的数组序号只能从1开始endc(k+1)=sqrt(a(k+1)^2+b(k+1)^2);endsubplot(211);plot(t,x);title('原始信号');xlabel('时间/s'); f=(0:m-1)/(N*dt);subplot(212);plot(f,c);title('Fourier变换’)；xlabel('频率/Hz');ylabel('振幅') （2）快速傅里叶变换fs=input('please input the fs:');%设定采样频率N=input('please input the N:');%设定数据长度(t,x)(先导入t,x的值)subplot(211);plot(t,x);%作正弦信号的时域波形axis([0,1,-1,1]);%设置坐标轴长度title('正弦信号时域波形');%进行FFT变换并做频谱图y=fft(x,N);%进行fft变换mag=abs(y);%求幅值f=(0:N-1)*fs/N;%横坐标频率的表达式为f=(0:M-1)*Fs/M; subplot(212);plot(f,mag);%做频谱图title('正弦信号幅频谱图');2.数据导入记事本数据导入代码方法一：x = load('E:\文件名.txt')plot(x(:, 1), x(:, 2))方法二：x=importdata(文件名)x1=x.data(行,列)%例如x1=x.data(1:100,2)(记事本中只有数据时，代码)x=importdata(文件名)x1=x(行,列)3.转速数据处理练习（1）转速处理代码离散变换a=importdata('speed1.txt')t=a(1:1000,1)x=a(1:1000,2)N=1000;dt=0.001;n=0:N-1; t=n*dt;%序号序列和时间序列m=floor(N/2)+1;%分解a，b的最大序号值为分解的N/2个参数加参数a0f=(0:m-1)/(N*dt);a=zeros(1,m);b=zeros(1,m);%产生a，b两个为零的序列for k=0:m-1,for i=0:N-1a(k+1)=a(k+1)+2/N*x(i+1)*cos(2*pi*k*i/N);b(k+1)=b(k+1)+2/N*x(i+1)*sin(2*pi*k*i/N);%MATLAB中的数组序号只能从1开始endc(k+1)=sqrt(a(k+1)^2+b(k+1)^2);endsubplot(211);plot(t,x);title('原始信号');xlabel('时间/s');subplot(212);plot(f,c);axis([0,20,0,2]);title('Fourier变换’)；xlabel('频率/Hz');ylabel('振幅')频率1，4,8,17HzFFTa=importdata('speed1.txt')t=a(1:1000,1)x=a(1:1000,2)N=1000;fs=1000subplot(211);plot(t,x);%进行FFT变换并做频谱图y=fft(x,N);%进行fft变换mag=abs(y);%求幅值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(212);plot(f,mag);%做频谱图title('正弦信号幅频谱图');axis([0,20,0,400]);频率1，4,8,17Hz改进clc;clear;a=importdata('speed1.txt')t=a(1:1000,1)x=a(1:1000,2)N=1000;fs=1000subplot(311);plot(t,x);title('原始信号')%进行FFT变换并做频谱图y=fft(x,N);%进行fft变换mag=abs(y);%fft模值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(312);plot(f,mag);%做频谱图title('fft模值—频率');axis([0,20,0,400]);yy=mag/(N/2);yy(1)=mag(1)/N;subplot(313);plot(f,yy);axis([0,20,0,2]);title('幅—频');角度形式clc;clear;a=importdata('speed1.txt') t=a(1:1000,1)x=a(1:1000,2)xd=x*360/60N=1000;fs=1000subplot(411);plot(t,x);title('原始信号') subplot(412);plot(t,xd);title('原始信号')%进行FFT变换并做频谱图y=fft(xd,N);%进行fft变换mag=abs(y);%fft模值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(413);plot(f,mag);%做频谱图title('fft模值—频率');axis([0,20,0,2000]);yy=mag/(N/2);yy(1)=mag(1)/N;subplot(414);plot(f,yy);axis([0,20,0,10]);title('幅—频');弧度形式clc;clear;a=importdata('speed1.txt') t=a(1:1000,1)x=a(1:1000,2)xr=x*2*pi/60N=1000;fs=1000subplot(411);plot(t,x);title('原始信号')subplot(412);plot(t,xr);title('原始信号')%进行FFT变换并做频谱图y=fft(xr,N);%进行fft变换mag=abs(y);%fft模值;%横坐标频率的表达式为f=(0:M-1)*Fs/M; f=(0:N-1)*fs/Nsubplot(413);plot(f,mag);%做频谱图title('fft模值—频率');axis([0,20,0,50]);yy=mag/(N/2);yy(1)=mag(1)/N;subplot(414);plot(f,yy);axis([0,20,0,0.1]);title('幅—频');。
中国剩余定理（CRT）and扩展中国剩余定理（EXCRT）
中国剩余定理（CRT）and扩展中国剩余定理（EXCRT）今有物不知其数，三三数之余⼆；五五数之余三；七七数之余⼆。
问物⼏何？语⽂⽔平不⾼，⼤概翻译⼀下：今天Rothen钓了⼏个妹⼦，3个3个的数会余下2个，5个5个的数会余下3个，七个七个的数会剩下2个好吧，这样好像更难理解了,懒得翻译了，反正⼤家都看得懂，⼤佬们就先做会吧，反正对于您这种神犇那肯定是秒切啊！。
中国剩余定理好⾼级的东西啊，吓得我赶紧来个BFS( Baidu First Search)发现中国剩余定理⼜叫孙⼦定理，原因是孙⼦发明的（这名字绝了）然后对于上⾯的题⽬的解法为：(注意，中国剩余定理的解法只能⽤于⼏个模数两两互质的情况)1.⾸先找出5和7的 mod 3等于1的公倍数(70)，然后找出3和7的mod 5 等于1的公倍数(21),还有就是3 和 5的mod 7等于1的公倍数(15)2.答案就是上⾯求出来的三个数的乘积分别乘以"对应的余数" mod 所有模数的乘积例如 (70*2+21*3+15*2)mod (7*3*5)得到23纳尼！这么神奇的吗？真的就是这么神奇哦～本蒟蒻查阅对于这个定理的证明为：把上⾯的问题转化为多个个⼦问题：假设存在⼀个数x1 % 3 =2("%"就是取余数),那么x1就能表⽰为3k+2的形式(k >= 0)假设存在⼀个数x2 % 5 =3("%"就是取余数),那么x1就能表⽰为5k+3的形式(k >= 0)假设存在⼀个数x3 % 7 =2("%"就是取余数),那么x1就能表⽰为7k+2的形式(k >= 0)那么考虑如果存在⼀个x1 + x2 + x3使得它们满⾜上⾯的所有条件呢？⾸先有个特别基础的公式：A % B = C, 那么(A + B*K) % B =C这个应该很好理解，你可以这么想：跟上⾯的类似,A可以表⽰为若⼲倍的B加上C，⽽A+B*k就可以表⽰为：若⼲倍的B + C +K倍的B ，⾃然的对于(A+B*k)%B也是等于C了啊！回到刚才的问题，怎么样才能使得x1+x2+x3满⾜条件呢？若要使得(x1+x2+x3)%3仍然余下2，那么x2,x3就⼀定要是三的倍数若要使得(x1+x2+x3)%5仍然余下3，那么x1,x3就⼀定要是五的倍数若要使得(x1+x2+x3)%7仍然余下2，那么x1,x2就⼀定要是七的倍数所以问题就转化为了求出x1,x2,x3哦吼，从上⾯我们看出了x1,x2,x3的⼏个性质：1.x1是5和7的公倍数，⽽且x1 % 3 = 22.x2是3和7的公倍数，⽽且x2 % 5 = 33.x3是3和5的公倍数，⽽且x3 % 7 = 2同时⼜引⼊⼀个数学公式：如果a%b=c,那么(a*k)%b=(a%b+a%b+a%b+...+a%b(k个a%b) ) mod b=c×k mod b然后问题就可以⽤开篇的⽅法解决了呀！⽽且找出来的是最⼩满⾜条件的数。
有关数论算法-中国剩余定理
d = gcd(a, n)，或者无解。
2017/8/8
19
University of Science and Technology of China
模运算和模线性方程
20
University of Science and Technology of China
初等数论概念
最大公约数性质：
2017/8/8
6
University of Science and Technology of China
初等数论概念
互质数：如果gcd(a,b)=1,则称a与b为互质数；
如果两个整数中每一个数都与一个整数p互为质数，则它们的积 与p互为质数，即：
唯一因子分解：
注：这种解法需要整数的素因子分解，而素因子分解是一个很难的 问题（NP问题）
2017/8/8
8
University of Science and Technology of China
最大公约数
欧几里得算法
Th.31.9（GCD递归定理）： 对任何非负整数a和正整数b，有gcd(a,b) = gcd(b, a mod b) ； * 可以通过证明gcd(a,b)与 gcd(b, a mod b)能相互整除来证明该定理！ P526
3
3
3
-2
-11
3
3
3 3
1
0 1
-2
1 0
gcd(a, b) gcd(b, a modb) (d, x, y) = ( a, 1, 0 ) gcd(a,0) a
数值分析MATLAB代码
% LU分解% 算例A=[1 2 3;1 3 5;1 3 6]，b=[2 3 4]'function lufun(A,b)n=length(A);L=zeros(n);U=zeros(n);y=ones(n,1);x=ones(n,1);for i=1:nL(i,i)=1;endfor i=1:nfor j=1:nif i>jU(i,j)=0;L(j,i)=0;endendendfor i=1:nfor j=i:nU(i,j)=A(i,j);endfor j=i+1:nL(j,i)=A(j,i)/U(i,i);endfor j=i+1:nfor k=i+1:nA(j,k)=A(j,k)-L(j,i)*U(i,k);endendendfor k=1:ns=0;for j=1:k-1s=s+L(k,j)*y(j);endy(k)=b(k)-s;endfor k=n:(-1):1s=0;for j=k+1:ns=s+U(k,j)*x(j);endx(k)=(y(k)-s)/U(k,k);endL,U,y,x % 顺序Gauss消去法解线性方程组function x=gauss(A,b)n=length(b);for k=1:(n-1)m=A(k+1:n,k)/A(k,k);A(k+1:n,k+1:n)=A(k+1:n,k+1:n)-m*A(k,k+1:n);b(k+1:n)=b(k+1:n)-m*b(k);A(k+1:n,k)=zeros(n-k,1);endx=zeros(n,1);x(n)=b(n)/A(n,n);for k=n-1:-1:1x(k)=(b(k)-A(k,k+1:n)*x(k+1:n))/A(k,k)end% Jacobi-迭代法% 系数矩阵，A=[5 2 1;-1 4 2;2 -3 10]% 常数向量，b=[-12;20;3]% 初值向量，x0=[-3;1;1]function Jacfun(A,b,x0)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=D\(L+U);F=D\b;format longn=1x=B*x0+Fwhile norm(x-x0)>=1e-3x0=x;n=n+1x=B*x0+Fif n>=500breakendend% SOR迭代法% A=[4 -1 0;-1 4 -1;0 -1 4];% b=[1;4;-3];% x0=[0 0 0]';function sorfun(A,b,x0)w=0.8;D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);n=1x=inv(D-w*L)*((1-w)*D+w*U)*x0+w*inv(D-w*L)*bwhile norm(x-x0)>=1e-6x0=x;n=n+1x=inv(D-w*L)*((1-w)*D+w*U)*x0+w*inv(D-w*L)*bif n>=500breakendend% GS-迭代法% 系数矩阵，A=[5 2 1;-1 4 2;2 -3 10]% 常数向量，b=[-12;20;3]% 初值向量，x0=[-3;1;1]function gsfun(A,b,x0)D=diag(diag(A));L=-tril(A,-1);U=-triu(A,1);B=(D-L)\U;F=(D-L)\b;n=1x=B*x0+Fwhile norm(x-x0)>=1e-9x0=x;n=n+1x=B*x0+Fif n>=500breakendend% 乘幂法求主特征值特征向量% v特征值,u特征向量% m----最大特征值% 算例A=[-12 3 3;3 1 -2;3 -2 7]function cmfun(A)k=length(A);n=0;m1=0;u=ones(k,1);while n<=500n=n+1v=A*ufor i=1:n[m,i]=max(abs(v));endm=v(i);u=v/mif abs(m-m1)<1e-8breakendm1=m;end% QR分解% 算例A=[1 2 2;2 2/3 -1;2 1/3 -1]，求A的QR分解function qrfun(A)n=length(A);A0=A;s=0;for j=1:n-1for i=j:ns=s+A(i,j)^2;enda=sign(A(j,j))*sqrt(s);b=A(j,j)+a;u=A(j:n,j);u(1)=b;t=length(u);R1=eye(t)-u*u'/(a*b);H=blkdiag(eye(n-t),R1);A1=H*A;A=A1;s=0;endR=AQ=A0/REnd 用雅可比法求特征值和特征向量算例A=[1 1 0.5;1 1 0.25;0.5 0.25 2];function jcbfun(A)n=length(A);Q0=eye(n);I=1;J=1;for k=1:400B=triu(A,1)+tril(A,-1);% 求矩阵B中绝对值最大的元素及所在的行和列a=0;for i=1:nfor j=1:nif a<abs(B(i,j))I=i;J=j;a=abs(B(i,j));endendendI;J;a=B(I,J);% 求平面旋转变换矩阵的相位角tif A(I,I)==A(J,J)t=pi*sign(A(I,J))/4;elsed=-2*A(I,J)/(A(I,I)-A(J,J));t=atan(d)/2;end% 求平面旋转矩阵PP=eye(n);P(I,I)=cos( t);P(I,J)=sin( t);P(J,I)=-sin( t);P(J,J)=cos( t);kA1=P'*A*P;A=A1Q=Q0*PQ0=Q;B=triu(A,1)+tril(A,-1);if norm(B,'fro')<=1e-9breakendend% A为实对称矩阵时，则可约化为三对角矩阵% 算例A=[1 2 1 2;2 2 -1 1;1 -1 1 1;2 1 1 1];function sdjfun(A)a1=sign(A(2,1))*sqrt(A(2,1)^2+A(3,1)^2+A(4,1)^2);b1=A(2,1)+a1;u1=[b1;A(3,1);A(4,1)];R1=eye(3)-u1*u1'/(a1*b1);H1=blkdiag(eye(1),R1)A1=H1*A*H1A=A1;a2=sign(A(3,2))*sqrt(A(3,2)^2+A(4,2)^2);b2=A(3,2)+a2;u2=[b2;A(4,2)];R2=eye(2)-u2*u2'/(a2*b2);H2=blkdiag(eye(2),R2)A2=H2*A*H2H=H2*H1A=A2% 拉格朗日插值函数% 算例x=[-2 -0.8 0.4 1.2];y=[1 1.4 1.8 2.0];function lagfun(x,y)syms t;n=length(x);s=0;for i=1:nla=y(i);for j=1:nif j~=ila=la*(t-x(j))/(x(i)-x(j));endends=s+la;simplify(s);ends=subs(s,'t','x');s=collect(s);s=vpa(s,6)% 牛顿插值法% 算例x=[-1 0 1 2 3];y=[-2 1 3 4 8]; function newton(x,y)syms p;s=y(1);t=0;d=1;n=length(x);for i=1:n-1for j=i+1:nt(j)=(y(j)-y(i))/(x(j)-x(i));endtemp(i)=t(i+1);d=d*(p-x(i));s=s+temp(i)*d;y=t;ends=subs(s,'p','x');y=simplify(s)龙贝格求积function romberg(a,b,f)tol=0.5e-5;n=1;h=b-a;delt=1;x=a;k=0;R=zeros(4,4);R(1,1)=h*(f(a)+f(b))/2;while delt>tolk=k+1h=h/2;s=0;for j=1:nx=a+h*(2*j-1);s=s+f(x);endR(k+1,1)=R(k,1)/2+h*s;n=2*n;for i=1:kR(k+1,i+1)=((4*i)*R(k+1,i)-R(k,i))/(4*i-1);enddelt=abs(R(k+1,k)-R(k+1,k+1)); s=R(k+1,k+1)endfunction z=f(x)z=@(x)4/(1+x^2); % 三次样条插值（第一类边界条件）% X为插值点，Y为插值点对应的函数值，dY为插值点两端点对应的一阶导数值% 算例X=[1 2 3];Y=[2 4 2];dY=[1 -1];function spline3(X,Y,dY)N=size(X,2);m=5;s0=dY(1); sN=dY(2);h=zeros(1,N-1);for i=1:N-1h(1,i)=X(i+1)-X(i);endd(1,1)=6*((Y(1,2)-Y(1,1))/h(1,1)-s0)/h(1,1);d(N,1)=6*(sN-(Y(1,N)-Y(1,N-1))/h(1,N-1))/h(1,N-1);for i=2:N-1d(i,1)=6*((Y(1,i+1)-Y(1,i))/h(1,i)-(Y(1,i)-Y(1,i-1))/h(1,i-1))/(h(1,i)+h(1,i-1));endmu=zeros(1,N-1);md=zeros(1,N-1);md(1,N-1)=1;mu(1,1)=1;for i=1:N-2u=h(1,i+1)/(h(1,i)+h(1,i+1));mu(1,i+1)=u;md(1,i)=1-u;endp(1,1)=2;q(1,1)=mu(1,1)/2;for i=2:N-1p(1,i)=2-md(1,i-1)*q(1,i-1);q(1,i)=mu(1,i)/p(1,i);endp(1,N)=2-md(1,N-1)*q(1,N-1);y=zeros(1,N);y(1,1)=d(1)/2;for i=2:Ny(1,i)=(d(i)-md(1,i-1)*y(1,i-1))/p(1,i);endx=zeros(1,N);x(1,N)=y(1,N);for i=N-1:-1:1x(1,i)=y(1,i)-q(1,i)*x(1,i+1);endM=x;syms t;digits (m);for i=1:N-1pp(i)=M(i)*(X(i+1)-t)^3/(6*h(i))+M(i+1)*(t-X(i))^3/(6*h(i))+(Y(i)-M(i)*h(i)^2/6)*(X(i+1)-t)/h(i)+(Y(i+1)-M(i+1)*h(i)^2/6)*(t-X(i))/h(i)pp(i)=simplify(pp(i));f=sym2poly(pp(i));if length(f)~=4tt=f(1:3); f(1:4)=0; f(2:4)=tt;ends=sym(f,'d');s=poly2sym(s,'t');fprintf('在区间[%f,%f]内',X(i),X(i+1));send% 四级四阶龙格库塔方法% 算例y'=1-y ,y(0)=0（0<=x<=1）,步长h=0.1% Runge_Kutta( 0, 0, 0.1, 20)function Runge_Kutta(x0, y0, h, n)x=zeros(n+1,1);y=zeros(n+1,1);x(1)=x0;y(1)=y0;for i=1:nx(i+1)=x(i)+h;k1=h*f(x(i),y(i));k2=h*f(x(i)+1/2*h,y(i)+1/2*k1);k3=h*f(x(i)+1/2*h,y(i)+1/2*k2);k4=h*f(x(i)+h,y(i)+k3);y(i+1)=y(i)+1/6*(k1+2*k2+2*k3+k4);endT=[x,y]% 显式欧拉法function E=Euler(x0,y0,xn,n)x=zeros(n+1,1);y=zeros(n+1,1);x(1)=x0;y(1)=y0;h=(xn-x0)/n;for i=1:nx(i+1)=x(i)+h;y(i+1)=y(i)+h*f(x(i),y(i));endT=[x,y]。