Register to reply

MATLAB code has stopped working

by hunt_mat
Tags: code, matlab, stopped, working
Share this thread:
hunt_mat
#1
Jan21-13, 09:08 AM
HW Helper
P: 1,583
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?
Phys.Org News Partner Science news on Phys.org
Physical constant is constant even in strong gravitational fields
Montreal VR headset team turns to crowdfunding for Totem
Researchers study vital 'on/off switches' that control when bacteria turn deadly
Pythagorean
#2
Jan21-13, 07:33 PM
PF Gold
Pythagorean's Avatar
P: 4,293
Are you running it on the same computer?
hunt_mat
#3
Jan22-13, 05:06 AM
HW Helper
P: 1,583
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.

NemoReally
#4
Jan22-13, 07:43 AM
P: 194
MATLAB code has stopped working

Quote Quote by hunt_mat View Post
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.
Can you get somebody else to run it or post it here for somebody to run / look at?
hunt_mat
#5
Jan22-13, 08:06 AM
HW Helper
P: 1,583
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


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