Dismiss Notice
Join Physics Forums Today!
The friendliest, high quality science and math community on the planet! Everyone who loves science is here!

MATLAB code has stopped working

  1. Jan 21, 2013 #1

    hunt_mat

    User Avatar
    Homework Helper

    A few months ago my hard drive crashed and I managed to recover the data and I tried to run some of my old programs (one solving an ODE via Newton's method. Before the crash, the error was 10^-8 after around 3 iterations, now it won't converge. It looks fine, ithe than writing the code from scratch again, I have no idea what to do.

    Suggestions?
     
  2. jcsd
  3. Jan 21, 2013 #2

    Pythagorean

    User Avatar
    Gold Member

    Are you running it on the same computer?
     
  4. Jan 22, 2013 #3

    hunt_mat

    User Avatar
    Homework Helper

    I am.

    The program had to be recovered from the failed hard drive, I am thinking that could be the reason why. It could be corrupted.
     
  5. Jan 22, 2013 #4
    Can you get somebody else to run it or post it here for somebody to run / look at?
     
  6. Jan 22, 2013 #5

    hunt_mat

    User Avatar
    Homework Helper

    I tried to run it on another computer, previously it ran and converged within 4 iteration to an error of 10^-8, now it won't converge for some reason.

    The code:
    clear;
    L=201; %Defines the number of steps, must always be an integer
    x=linspace(-6,6,L); %Defines x
    B=0.2; %This is the Bond number.
    p=0.25*B*exp(-x.^2); %The pressure distribution
    E_b=0.2; %This is the electric Bond number.
    f=zeros(1,L+1);
    f(1,1:L)=0.01*exp(-x.^2); %have some small initial guess for the wave profile
    choice=false; %This selects the option to fix the Froude number of the minimum of the wave.
    Q=-0.6; %This is the Froude number (choice=true), or the minimum of the profile (choice=false)
    f(L+1)=Q; %this is the Foude number, either an exact value or an initial guess.
    %The way forward is to use Newton's method to reduce the equation down to a linear algebra problem and then use Octave
    %to solve the linear algrbra problem.

    %The first step is to compute the Jacobian matrix
    %Do this by diferencing between to close points and dividing through by that difference
    df=10^-10;
    J=zeros(L+1,L+1); %This will be the Jacobian.
    X=zeros(L+1,1);
    error=0.5;
    while (error>10^-10)
    for i=1:L+1
    for j=1:L+1
    T=zeros(1,L+1); %Set all the entries to zero again
    T(j)=df;
    J(i,j)=(g(x,f+T,p,Q,i,choice,B,E_b)-g(x,f-T,p,Q,i,choice,B,E_b))/(2*df); %Calculates the (i,j)th element of the Jacobian
    end
    end
    for m=1:L+1
    X(m)=-g(x,f,p,Q,m,choice,B,E_b);
    end
    %Now solve the system of equations
    sol=J\X;
    error=norm(X)
    f=f+sol';
    end

    %Now plot the solution!
    plot(x,f(1:L));
    xlabel('x');
    ylabel('f(x-ct)');

    Using:

    function r=g(x,f,p,Q,i,flag,B,E_b)
    N=length(x);
    h=1;
    a_1=1-f(N+1);
    a_2=h^4/90;
    a_3=1.5/h;
    a_4=-0.5*(B-1/3)*h^2;
    a_5=0.5*E_b*h;
    dx=x(2)-x(1);
    ff=f(1:N);
    f_ghost=[0 0 ff 0 0]; %Add in the ghost cells for f
    %Prepare for the Hilbert transform
    x_half=zeros(1,N-1);
    f_half=zeros(1,N-1);
    for k=1:N-1
    x_half(k)=0.5*(x(k)+x(k+1));
    f_half(k)=0.5*(ff(k)+ff(k+1));
    end
    grad=gradient(f_half,dx);
    if (i<N+1)
    hilbert=trapz(x_half,grad./(x_half-x(i)));
    end
    %Now comes the formulae for the finite differencing part of the algorithm.
    if (i<N+1)
    b_1=a_1*f_ghost(i+2);
    b_2=(a_2/dx^4)*(f_ghost(i+4)-4*f_ghost(i+3)+6*f_ghost(i+2)-4*f_ghost(i+1)+f_ghost(i)); %Fourth order derivative
    b_3=a_3*f_ghost(i+2)^2; %nonlinear term
    b_4=(a_4/dx^2)*(f_ghost(i+3)-2*f_ghost(i+2)+f_ghost(i+1)); %Second order derivative
    b_5=a_5*hilbert; %Hilbert transform term
    b_6=p(i); %Pressure term
    r=b_1+b_2+b_3+b_4+b_5+b_6;
    else
    if (flag==true)
    r=f(N+1)-Q;
    else
    r=f(1+0.5*(N-1))-Q;
    end
    end
     
Know someone interested in this topic? Share this thread via Reddit, Google+, Twitter, or Facebook




Similar Discussions: MATLAB code has stopped working
  1. Matlab codes (Replies: 2)

  2. Matlab code (Replies: 0)

  3. Error code in Matlab (Replies: 3)

Loading...