matlab一个很小的编程题目,求助!!!

本人matlab基本上是不大懂的。。。本来想用C语言实现,结果发现C语言在数学问题上实在有限,题目不难,但是实在不会用matlab把代码写出来
因此求助各位大神,耽误您一点小小的时间帮我写一下这个题目的matlab代码行吗?谢谢各位!我会附上我的c语言代码,如果有需要,可以直接用。

题目如下:一个小蚂蚁在平面格点上移动,1/4概率上下左右,每次移动一步,比如从(0,0)到(0,1),不能斜着移动。
设总共走了n步,
记这n步中,走到过的格点数目为R(重复走到的格点,只能算一个)
记这n步中,到过int(ln(n))次的格点数目为r
计算当n趋于无穷时 r/R的比值

例如:走了n=10步,可能只到过8不同的格点(有重复走到的格点),此时R=8
在n=10步中,int(ln(10))=2,此时统计到过两次的格点数目,可能有1个,r=1

我的c语言代码如下
#include<cmath>
#include <ctime>
#include <cstdlib>
using namespace std;
long const N=100;
long const M=2*N;
int main()
{
int a[M+1][M+1]={0};
int b;
int j=N,k=N;
double R=0,r=0;
srand(unsigned(time(0)));
for(int i=1;i<=N;i++)
{
b=1+4*rand()/RAND_MAX;
switch (b)
{
case 1:j++;a[j][k]=a[j][k]+1;break;
case 2:k++;a[j][k]=a[j][k]+1;break;
case 3:j--;a[j][k]=a[j][k]+1;break;
case 4:k--;a[j][k]=a[j][k]+1;break;
default:break;
}
}
for(int i=0;i<=M;i++)
for(int s=0;s<=M;s++)
{
if(a[i][s]!=0)
R=R+1;
}
if(int(log(double(N)))!=0)
{
for(int i=0;i<=M;i++)
for(int s=0;s<=M;s++)
{
if(a[i][s]==int(log(double(N))))
r=r+1;
}
}
cout<<r<<endl<<R<<endl;
cout<<r/R<<endl;
return 0;
}

楼上两位的回答都很用心,也很精彩,赞一个。

 

我的代码主要有以下优点:

(1)用稀疏矩阵存储a,克服内存不足问题(N取100万,使用的内存还不到20M)。

(2)绘图动态显示N次模拟过程中r/R的变化。

 

代码如下(同时已作为附件上传):

N = 1000000;
M = 2*N;
a = sparse(M+1, M+1);
j = N + 1;
k = N + 1;
b = ceil(4*rand(1,N));

% 绘图显示计算过程(为提高效率,每n次循环输出一个点)
n = 500;
v = zeros(1, fix(N/n)+1) * NaN;
h = plot(1 : fix(N/n)+1, v, 'b-', NaN, NaN, 'ro');
xlabel('N');
ylabel('r/R');
set(gcf, 'DoubleBuffer', 'on');
set(gca, 'Xlim', [1 fix(N/n)+1]);
t=ceil(exp(1:15));
set(gca,'xtick',t/n,'xgrid','on','xticklabel',t)

for i = 1:N
    switch b(i)
        case 1, j=j+1; a(j,k)=a(j,k)+1;
        case 2, k=k+1; a(j,k)=a(j,k)+1;
        case 3, j=j-1; a(j,k)=a(j,k)+1;
        case 4, k=k-1; a(j,k)=a(j,k)+1;
    end
    % 更新绘图
    if ~rem(i, n) || i == N
        % 注意:不能用 sum(a(:)~=0) 进行计算,否则容易导致内存不足
        R = full(sum(sum(a~=0)));
        c = fix(log(i));
        r = full(sum(sum(a==c)));
        idx = fix(i/n) + 1;
        v(idx) = r/R;
        set(h(1), 'yData', v);
        set(h(2), 'xData', idx, 'yData', v(idx));
        title(['N = ' int2str(i)]);
        drawnow
    end
end

% 输出结果
fprintf('R=%i, r=%i, r/R=%.3g\n\n', R, r, r/R);
% 统计各方向移动的次数(验证随机数的均匀性)
for i = 1 : 4
    fprintf('方向%i的次数为%i\n', i, sum(b==i));
end

某次程序的输出如下:

R=204146, r=2856, r/R=0.014
方向1的次数为249327
方向2的次数为250179
方向3的次数为250085
方向4的次数为250409

 

简单说明几点:

1、由于是随机模拟,每次运行的结果都会有差别。

2、移动是一个前后关联的过程,所以随机数序列不仅要求均匀,还应该独立(相邻的随机数之间不相关)。

3、图中的虚线表示 int(ln(n)) 发生变化的N,计算格点次数变了,所以通常表现为不连续。

4、从图中的变化趋势看,没有收敛到某一个数的明显迹象,我运行两次的结果分别是0.0118和0.014。

5、在我的机器上,取N=100万,运行一次的时间大约是10分钟。

6、程序对MATLAB版本没有特别要求,在6.5、2008a、2012b上测试过,都可以正常运行。

追问

看了你的结果,发现这个题被我想简单了,我原以为极限存在,那么接着就好做了,但既然不存在,我就又不知道该咋办了。本来题目是这样的:记这n步中,走到过的格点数为R(重复到的格点,算一个);记这n步中,到过int(x*ln(n))次的格点数为r(x);计算当n趋于无穷时r(x)/R的比值,r(x)/R是x的函数,画出r(x)/R的图像(关于x)能把代码和图像都弄出来吗?如果能把这个帮我解决了,分就给你

追答

你没给出x的定义域,我根据情况试着取了1-4之间。N不可能无限大,我取了20万,绘图仍然是动态更新的,但到最后只保留最终的一组曲线如下。从曲线看有点像指数曲线,但不确定。

代码如下:

N = 200000;
M = 2*N;
a = sparse(M+1, M+1);
j = N + 1;
k = N + 1;
b = ceil(4*rand(1,N));
 
% 绘图显示计算过程(为提高效率,每50次循环输出一个点)
n = 500;
v = zeros(1, fix(N/n)+1) * NaN;
x = 1:.1:4;
rx = zeros(size(x)) * NaN;
h = plot(x,rx);
xlabel('x');
ylabel('r(x)/R');
set(gcf, 'DoubleBuffer', 'on');
 
for i = 1:N
    switch b(i)
        case 1, j=j+1; a(j,k)=a(j,k)+1;
        case 2, k=k+1; a(j,k)=a(j,k)+1;
        case 3, j=j-1; a(j,k)=a(j,k)+1;
        case 4, k=k-1; a(j,k)=a(j,k)+1;
    end
    % 更新绘图
    if ~rem(i, n) || i == N
        % 注意:不能用 sum(a(:)~=0) 进行计算,否则容易导致内存不足
        R = full(sum(sum(a~=0)));
        c = log(i);
        for ii=1:length(x)
            rx(ii) = full(sum(sum( a==fix(x(ii)*c) )));
        end
        set(h, 'yData', rx/R);
        title(['N = ' int2str(i)]);
        drawnow
    end
end

温馨提示:答案为网友推荐,仅供参考
第1个回答  2013-07-22
一、你的代码已经非常清晰,只不过存在不同编程方法的转换问题。
PS:①关于随机数生成问题:matlab中的对应函数是rand,产生0~1之内的均匀分布的数;② 向下取整函数:C语言中的int函数用Matlab中的floor函数替换。
二、M文件代码:
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%% 2013-07-22:改写C代码,edited by ljwgbb
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear %%%% 清屏及清除变量
%%%%%%% 1、参数设置
NumN=100;
NumM=NumN*2;
a=zeros(NumM+1,NumM+1);
b=0;
j=NumN;
k=NumN;
R=0;
r=0;
rand('state',20130722); %%%%% 随机数种子
%%%%%%% 2、算法部分
for i=1:1:NumN
%%%%% 产生满足均匀分布1~4的整数
temp=randperm(4);
b=temp(1,1);
switch b
case 1
j=j+1;a(j,k)=a(j,k)+1;
case 2
k=k+1;a(j,k)=a(j,k)+1;
case 3
j=j-1;a(j,k)=a(j,k)+1;
case 4
k=k-1;a(j,k)=a(j,k)+1;
otherwise
warning('Unexpected Error:b is not interage.');
end
end

for i=1:1:NumM+1
for s=1:1:NumM+1
if a(i,s)>0
R=R+1;
end
end
end
MyVal=floor(log(NumN))
if MyVal~=0
for i=1:1:NumM+1
for s=1:1:NumM+1
if a(i,s)==MyVal
r=r+1;
end
end
end
end
%%%%%%% 3、结果部分
MyRes=[r,R,r/R]

最终的运行结果如下:
MyRes =

4.0000 46.0000 0.0870

三、编写成函数
function MyRes=MyFunV1(NumN)
%% 数值模拟函数
NumM=NumN*2;
a=zeros(NumM+1,NumM+1);
b=0;
j=NumN;
k=NumN;
R=0;
r=0;
rand('state',20130722); %%%%% 随机数种子
%%%%%%% 2、算法部分
for i=1:1:NumN
%%%%% 产生满足均匀分布1~4的整数
temp=randperm(4);
b=temp(1,1);
switch b
case 1
j=j+1;a(j,k)=a(j,k)+1;
case 2
k=k+1;a(j,k)=a(j,k)+1;
case 3
j=j-1;a(j,k)=a(j,k)+1;
case 4
k=k-1;a(j,k)=a(j,k)+1;
otherwise
warning('Unexpected Error:b is not interage.');
end
end

for i=1:1:NumM+1
for s=1:1:NumM+1
if a(i,s)>0
R=R+1;
end
end
end
MyVal=floor(log(NumN))
if MyVal~=0
for i=1:1:NumM+1
for s=1:1:NumM+1
if a(i,s)==MyVal
r=r+1;
end
end
end
end
%%%%%%% 3、结果部分
MyRes=[r,R,r/R]追问

那你这个程序的运行结果是什么呀,我的意思是,题目是要求计算n趋于无穷的时候,那个比值是多少,那这个结果是多少啊?是0.0870吗。。。?由于本人matlab啥也不懂,所以还有个疑问,就是你的“二、M文件代码”和“三、编写成函数”这俩有啥区别,我感觉怎么差不多呢,为啥你列了俩。。。?

追答

一、M文件,相当于语句集合,其中参数是确定的,如NumN=200;函数,相当于功能模块,可对输入参数有效控制。现使用函数,实现你的追问:

二、M代码如下:

%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%% 调用MyFunV1求解

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clc,clear

NumN=[100,500,1000:1000:8000];

NumOpe=numel(NumN);

MyRes=ones(NumOpe,6);     %%%%%[ i1,NumN,r ,R,r/R,RunTime]

for i1=1:NumOpe

    i1

    tic

    MyRes0 = MyFunV1( NumN(1,i1) );

    RunTime=toc;

    MyRes(i1,:)=[i1,NumN(1,i1),MyRes0,RunTime];

end

disp(sprintf('\nCongratulations!'))

三、运算结果如下:

最终近似结果为:0.0227,但趋势并不明显。此时数组不能再进一步扩充!需改用其它算法或其它存储模式。


四、分析:上述涉及到数组大小容量的问题,不能无限变大,这也许就是这道题的难点和重点吧:

关于Matlab的矩阵大小可参考命令:memory

算法需要改变,将NumM×NumM的矩阵规模降为NumM×K(K<5):popo0104提供的C++算法更为合理(貌似改编的Matlab代码有些小问题)。


本回答被网友采纳
第2个回答  2013-07-22

C++版本:

#include <iostream>
#include <cmath>
#include <ctime>
#include <cstdlib>
using namespace std;

void SetB(long i, int a[][2], int b[][3], int *n);

int main(int argc, char *argv[])
{
const long N = 50000;  //总次数 
double r=0.0, R=0.0;

//int a[N][2] = {0}, b[N][3]={0};
int (*a)[2] = new int[N][2];  //按次序记录点的坐标 
int (*b)[3] = new int[N][3];  //记录走过的点,前两列是坐标,第三列是经过该点的次数 
int t, s, np;

a[0][0] = 0;a[0][1] = 0;      //假设初始位置(0,0) 
b[0][0] = 0;a[0][1] = 0;b[0][2] = 1; //经过的第一个点(0,0),记录为1次 
np = 1;                              //经过的不重复的点数目 

srand(unsigned(time(0)));
for(long i=1; i<N; i++)
{
//cout << rand() << "  " << RAND_MAX <<endl;
t=1+4*rand()/(RAND_MAX + 1.0); //产生均匀随机数 1,2,3,4 
//cout << "t=" << t << endl;
switch (t)
{
case 1://向左移 
{a[i][0]=a[i-1][0] - 1; a[i][1]=a[i-1][1]    ; SetB(i,a, b, &np); break;}
case 2://向上移 
{a[i][0]=a[i-1][0]    ; a[i][1]=a[i-1][1] + 1; SetB(i,a, b, &np); break;}
case 3://向右移 
{a[i][0]=a[i-1][0] + 1; a[i][1]=a[i-1][1]    ; SetB(i,a, b, &np); break;}
case 4://向下移 
{a[i][0]=a[i-1][0]    ; a[i][1]=a[i-1][1] - 1; SetB(i,a, b, &np); break;}
default:cout << t << "error!\n";
}
if(i%10000 ==0 ) cout << "Now i = "<< i << ", np= " << np << endl;
}
/*for(long i=0;i<N;i++){
cout << "a[" << i << "][0]= " << a[i][0] <<" a[" << i << "][1]= " << a[i][1] << endl;
}
cout << "np = " << np << endl;
for(long i=0;i<np;i++){
cout << "b[" << i << "][0]= " << b[i][0] 
     <<" b[" << i << "][1]= " << b[i][1] 
 <<" b[" << i << "][2]= " << b[i][2]<< endl;
}*/
s = int(log(double(N)));
for(long i=0; i<np; i++){ //找出经过次数为 s 的点数目 
if(b[i][2] == s) r = r + 1;
}
R = (double)np;
cout << "N = " << N <<", int(log(double(N))) = " << s <<endl;;
cout << "r = " << r
     << ", R = " << R
 << ", r/R=" << r/R << endl;
 
delete [] a;
delete [] b;
return 0;
}
/*
每移动到一个位置,需要如下判断:
该点是否是已经经过的点,也就是在b中寻找是否有一样的点,如果有,那么经过次数加1,也就是b[j][2] + 1 ;
如果到达的点没有经历过,那么添加到b中,并且经历次数设置为1 
*/ 
void SetB(long k, int a[][2], int b[][3], int *n)
{
int sign = 0;//标记,看要检查的点是否经历过 
for(long j=0; j<k; j++){
if(a[k][0] == b[j][0] && a[k][1] == b[j][1])//如果经历过,则经历次数+1 
{
b[j][2] = b[j][2] + 1;
sign = 1;
break;
}
}
if(sign == 0){//这个点之前没有经历过,作为新点加入b中 
b[*n][0] = a[k][0];
b[*n][1] = a[k][1];
b[*n][2] = 1;
*n = *n + 1;
}

}

matlab:

%

clear;clc;

N=10000;

p=zeros(N,1);q = zeros(N,2);

t=zeros(N,1);

p(1) = 0 + 0*i;q(1,1) = p(1);q(1,2) = 1;

npoint = 1; %

for k=2:N

    t(k)=rand();

    if(t(k) < 1/4)     %左移

        p(k) = p(k-1) - 1;

        [npoint, q] = SetQ(k,p,q,npoint);

    elseif(t(k) < 1/2) %上移

        p(k) = p(k-1) + i;

        [npoint, q] = SetQ(k,p,q,npoint);

    elseif(t(k) < 3/4) %右移

        p(k) = p(k-1) + 1;

        [npoint, q] = SetQ(k,p,q,npoint);

    else               %下移

        p(k) = p(k-1) - i;

        [npoint, q] = SetQ(k,p,q,npoint);

    end

end

R = npoint;

n1 = floor(log(N))

r = 0;

for k = 1:npoint

    if (q(k,2) == n1)

        r = r + 1;

    end

end

r,R

r_R = r/R


%---------------------------------------------------

%SetQ.m

function [npoint, q] = SetQ(k, p, q, npoint)


sign = 0;

for j = 1:npoint  %判断是否已经走过点,若是,统计加1

    if(p(k) == q(j,1))

        q(j,2) = q(j,2) + 1;

        sign = 1;

    end

end

if(sign == 0)    %若没有走过,则添加入新点

    npoint = npoint + 1;

    q(npoint,1) = p(k);

    q(npoint,2) = 1;

end


        



追问

你这个有运行结果吗?我的意思题目是要求计算n趋于无穷的时候,那个比值是多少,那你的这个程序运行的结果是多少?

追答

我写的这个能算到很大,关键是你内存足够

第3个回答  2013-07-22
好巧,我研究酒鬼回家这个数学问题的时候写过差不多的程序,可惜代码找不到了,而且不是matlab,没法发给你!不过按照你的程序思路,貌似用不用matlab都一样啊?而且matlab写没有c语言直观
相似回答