Thread Subject:

Matrix is close to singular or badly scaled

Hi all,

I have a problem in solving a pde function ur*(partial delta/partial x)^2-ux*(partial^2) delta/(partial ^2)x+2*sqrt(ur*P)=b*sin(phi-delta).
The value of phi=0.5*exp(x-Ng/2)^2.(Ng=64)
I use the central difference for partial delta/ partial x=(delta(i+1)-delta(i-1))/(2*h) and (partial^2) delta/(partial ^2)x=delta(i+1)-2*delta(i)+delta(i-1)/h^2.

But matlab says Matrix is close to singular or badly scaled and cannot give the answer. So I modify the equation by add (epsilon/t+b)*delta(i) at both sides, then I get the figure and if I change the value of epsilon to 50, the amplitude of the figure decreases, just like epsilon is a damping factor.

My questions are as follows. Why by adding something at both sides, the equation can be solved? Theorectically the change of epsilon won't affact the result since it appears both sides. Why it plays a role of a damping factor?

Thank you!!
my codes is attached.
myfunction.m
global Ng ux ur b P m t h epsilon x N H
Ng=64;
ux=6;
ur=1;
b=1*ux;
P=0.01;
t=0.1; % Time Steps;
h=0.1; % Spacial Spacing;
epsilon=0.5;
x=linspace(0,Ng-h,Ng/h)';
N=length(x);
A=zeros(N);
for counter=2:N-1; % Center for Delta
   A(counter,counter)=2*ux/h^2+epsilon/t+b;
   A(counter,counter-1)=-ux/h^2;
   A(counter,counter+1)=-ux/h^2;
end;
A(1,1)=2*ux/h^2+epsilon/t+b;
A(1,2)=-ux/h^2;
A(1,N)=-ux/h^2;


A(N,N)=2*ux/h^2+epsilon/t+b;
A(N,N-1)=-ux/h^2;
A(N,1)=-ux/h^2;

H=inv(A);
phi0=0.5*exp(-0.1*(x-(Ng)/2).^2); % Disturbance in Phi
delta_guess=zeros(N,1);
delta0=algpde(phi0,delta_guess); % Corresponding Delta
plot(x,delta0)
alg.m
function delta=algpde(phi,delta_guess);
% PDE model for 64 generator ring system with internal impedance

global Ng ux ur b P m t h epsilon x N H A

error=[];
delta_old=delta_guess;
for i=1:25;
for counter=2:N-1;
   y(counter,1)=(epsilon/t+b)*delta_old(counter)-ur/(4*h^2)*(delta_old(counter+1)-delta_old(counter-1))^2+...
                b*sin(phi(counter)-delta_old(counter))-...
                2*sqrt(ur*P)/(2*h)*(delta_old(counter+1)-delta_old(counter-1));
end;
y(1,1)=(epsilon/t+b)*delta_old(1)-ur/(4*h^2)*(delta_old(2)-delta_old(N))^2+...
       b*sin(phi(1)-delta_old(1))-...
       2*sqrt(ur*P)/(2*h)*(delta_old(2)-delta_old(N));
y(N,1)=(epsilon/t+b)*delta_old(N)-ur/(4*h^2)*(delta_old(1)-delta_old(N-1))^2+...
       b*sin(phi(N)-delta_old(N))-...
       2*sqrt(ur*P)/(2*h)*(delta_old(1)-delta_old(N-1));
   
delta_new=H*y;
error=[error norm(delta_new-delta_old)];
delta_old=delta_new;
end;
delta=delta_new;
error

On 6/1/2012 6:01 PM, Shuang wrote:
> Hi all,
>
> I have a problem in solving a pde function ur*(partial delta/partial x)^2-ux*(partial^2)
>delta/(partial ^2)x+2*sqrt(ur*P)=b*sin(phi-delta).
> The value of phi=0.5*exp(x-Ng/2)^2.(Ng=64)
> I use the central difference for partial delta/ partial x=(delta(i+1)-delta(i-1))/(2*h)
>and (partial^2) delta/(partial ^2)x=delta(i+1)-2*delta(i)+delta(i-1)/h^2.
>
> But matlab says Matrix is close to singular or badly scaled and cannot give the answer.

What you need to do is first write down the full description of the pde in
normal math notation that one can read. The above is all mumble jumple to me
that is hard to read.

Also write down the domain and the initial/boundary conditions.

You might also want to check if the problem is well formed. But if
the teacher gave you thins, then I assume it is well formed already.

Then write down the finite difference scheme you used, and then write down what
the form of your A matrix that results when you discretize the pde at
each point, taking into account the boundary conditions
and how you discretized those. Boundary conditions are very important.

You can normally see the problem when you look at the A matrix form
when you do that. You do not need a computer for all of this. Paper
and pencil is all what you need.

That is what I had to do for my school HW's. We had to do all these
steps before even touching Matlab and coding anything. The coding
is the easy part. Do not jump to Matlab before you spend lots of
time in the analysis part of the problem and formulating the A
matrix and understating its form well.

If you do all that, and still have problem, then show all the above and
someone can help you more.

Many students seem to approach solving computation problems by keep
modifying code blindly until something 'works' and some numbers come out.

> My questions are as follows. Why by adding something at both sides, the equation
>can be solved? Theorectically the change of epsilon won't affact the result
>since it appears both sides. Why it plays a role of a damping factor?
>
> Thank you!!
> my codes is attached.
> myfunction.m

<snip code>

--Nasser

As far as I can see solving your equation without the epsilon-term
means that the solution delta
is not unique: If delta is a solution, then so is delta + 2*pi*n for n
an arbitray integer.
Or which boundary conditions did you impose ?

Best wishes
Torsten.