1 编写函数:计算n的阶乘
% fact.m
function output = fact(n)
% FACT Calculate factorial of a given positive integer.
output = 1;
for i = 1:n
output = output*i;
end
2 基于蒙特卡洛法计算pi
%N=10000; 定义计算pi的函数并设置函数参数:重复次数N
function output=getpi(N)
n=0;
for i=1:N
x=rand(1);
y=rand(1);
if (x-0.5)^2+(y-0.5)^2<=0.5^2
n=n+1;
end
end
%fprintf('Pi = %d .\n',4*n/N);
output=4*n./N;
end
power=1:0.5:6;
pivalue=[];pierror=[];
NMC=30;
for i=1:length(power)
dummy=[];
for j=1:NMC
dummy=[dummy,getpi(10^power(i))];
end
pivalue=[pivalue,mean(dummy)];
pierror=[pierror,std(dummy)];
end
% plot(power,pivalue)
e=sqrt(4*pi*(1-pi/4)./(10.^power));
figure(1);
errorbar(power,pivalue,e);
figure(2);
errorbar(power,pivalue,pierror);
% 结论:随着N的增大,误差逐渐变小,但边际效用递减。
3 解方程 log(x)=sin(x)
% solve log(x)=sin(x)
function output=solvelogxsinx(epsilonx)
epsilony=0.001;
x=2+epsilonx;
while(abs(log(x)-sin(x))>=epsilony)
x=x+epsilonx;
end
fprintf("The answer of log(x)=sin(x) is: %d ,\n and sin(x)=log(x)=%d.\n",x,sin(x));
output=x;
end
%plotsolvelogxsinx
power=1:6
epsilonx=1/100./(10.^(power/2));
solveresult=solvelogxsinx(epsilonx);
plot(epsilonx,solveresult)
% 结论:随着偏差的减小,最终数值解逐渐逼近真实解。
4 生成三阶幻方
a=zeros(3);
e=zeros(1,8);
b=1:9;
d=perms(b);
for i = 1:size(d,1)
a=reshape(d(i,:),3,3);
for j=1:3
e(j)=sum(a(j,:));
end
for j=1:3
e(3+j)=sum(a(:,j));
end
e(7)=a(1,1)+a(2,2)+a(3,3);
e(8)=a(1,3)+a(2,2)+a(3,1);
if e(1)==15 && e(2)==15 && e(3)==15 && e(4)==15 && e(5)==15 &&e(6)==15 &&e(7)==15 &&e(8)==15
disp(a)
disp('')
end
end
5 斐波那契数列
%fibonicci.m
function fib=fibonacci(n)
% 生成斐波那契数列
if n==1
fib=[1];
elseif n==2
fib=[1 1];
elseif n<1
fprintf('请输入大于1的整数\n');
else
fib=[1 1];
i=1;
while i<=n-2
fib(i+2)=fib(i)+fib(i+1);
i=i+1;
end
end
6 计算内部收益率
%getiir.m
function iir=getiir(T,C,A0)
% T:期数 C:现金流 A0:起始资金
if length(T)~=length(C)
fprintf('数据维度不一致\n');
else
eps=1000;
for iiri=0:0.001:1
sum=0;
for i=1:length(T)
sum=sum+C(i)/((1+iiri)^T(i));
end
if abs(sum-A0)<=eps
eps=abs(sum-A0);
iir=iiri;
end
end
end
7 贷款利率计算
%mortgage.m: get mortgage ...
function payment=mortgage(way,N,r,A0)
if way==1
%等额本息
sum=0;
for i=1:12*N
sum=sum+1/((1+r/12)^i);
end
payment=ones(1,12*N).*A0/sum;
elseif way==0
%等额本金
payment=ones(1,12*N);
for i=1:12*N
payment(i)=A0/12/N+(A0-A0*(i-1)/12/N)*r/12;
end
else
fprintf('月供方式输入错误!\n');
end
8 模拟自由落体运动
%fallinari.m
g=9.8;m=60;k=1.983 % 考虑空气阻力
x=0;v=0;t=0;
dt=0.1;dv=1;
while(dv>=0.01)
dv=(g-k*v/m)*dt;
v0=v;
v=v+dv*dt;
x=x+0.5*(v+v0)*dt;
fprintf('t=%f,v=%f,x=%f \n',t,v,x);
t=t+dt;
end
9 重复生日问题
function p=samebday(n)
dummy=1
for i=1:n
dummy=dummy*(365-(i-1))/365;
end
p=1-dummy
end
10 会面问题
甲、乙两人相约7点到8点在某地会面,先到者等候另一人20分钟,过时就可离去,试求这两人能会面的概率
% getmeeting.m
function output=getmeeting(N)
n=0;
for i=1:N
x=rand(1)*60;
if abs(x-y)<=20
n=n+1;
end
end
%fprintf('Pi = %d .\n',4*n/N);
output=n/N;
end
11 财富分配模拟
t_end=365*(65-18); % 365天*(65-18)年
ppl=100;%100人
fortune=zeros(t_end+1,ppl);
%初始化,100人,每人100元
fortune(1,:)=100;
money=fortune(1,:);
figure(1);
%hold on;
for t=1:t_end %每天随机分配一次
for i=1:ppl %每个人都要支出
% if money(i)>0 % 注释本行表示可贷款
% 如果没钱了,就不给钱了;但是不妨碍别人给自己钱
money(i)=money(i)-1; %自己的钱少一元:支出
picki=ceil(rand(1)*ppl); %随便选个人把支出的一元给出去
money(picki)=money(picki)+1; %那个幸运的人得到了1元
%end
end
fortune(t+1,:)=money; %随机分配一次之后保留数据
if t==round(t/100)*100 %为了加快速度,逢100画图
bar(sort(money)); % 排序后画图
title(['Day',num2str(t)]); % 在图像上标出天数
pause(0.00001); %pause的参数单位为秒。延迟一下,让计算等待画图完成
end
end