
#1
Jan2113, 09:08 AM

HW Helper
P: 1,584

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? 



#3
Jan2213, 05:06 AM

HW Helper
P: 1,584

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. 



#4
Jan2213, 07:43 AM

P: 194

MATLAB code has stopped working 



#5
Jan2213, 08:06 AM

HW Helper
P: 1,584

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,fT,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(xct)'); Using: function r=g(x,f,p,Q,i,flag,B,E_b) N=length(x); h=1; a_1=1f(N+1); a_2=h^4/90; a_3=1.5/h; a_4=0.5*(B1/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,N1); f_half=zeros(1,N1); for k=1:N1 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_halfx(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*(N1))Q; end end 


Register to reply 
Related Discussions  
Please I need a Help with matlab code (Controlling the LEGO NXT Using MatLab)  Math & Science Software  5  
Why isn't my algorithm code working?  Engineering, Comp Sci, & Technology Homework  4  
Matlab R2011a's new feature: "portable C/C++ code directly from MATLAB"  Math & Science Software  1  
COM Surrogate Has Stopped Working  Computing & Technology  2  
latex stopped working!  Forum Feedback & Announcements  2 