王半仙的学习分享

% fact.m

function output = fact(n) 
% FACT Calculate factorial of a given positive integer. 
output = 1; 
for i = 1:n
    output = output*i; 
end

%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的增大，误差逐渐变小，但边际效用递减。

% 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)
% 结论：随着偏差的减小，最终数值解逐渐逼近真实解。

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

%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

%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

%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

%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

function p=samebday(n)
    dummy=1
    for i=1:n
        dummy=dummy*(365-(i-1))/365;
    end
    p=1-dummy
end

% 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

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