Hi all,
I have a problem in solving a pde function ur*(partial delta/partial x)^2ux*(partial^2) delta/(partial ^2)x+2*sqrt(ur*P)=b*sin(phidelta).
The value of phi=0.5*exp(xNg/2)^2.(Ng=64)
I use the central difference for partial delta/ partial x=(delta(i+1)delta(i1))/(2*h) and (partial^2) delta/(partial ^2)x=delta(i+1)2*delta(i)+delta(i1)/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,Ngh,Ng/h)';
N=length(x);
A=zeros(N);
for counter=2:N1; % Center for Delta
A(counter,counter)=2*ux/h^2+epsilon/t+b;
A(counter,counter1)=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,N1)=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:N1;
y(counter,1)=(epsilon/t+b)*delta_old(counter)ur/(4*h^2)*(delta_old(counter+1)delta_old(counter1))^2+...
b*sin(phi(counter)delta_old(counter))...
2*sqrt(ur*P)/(2*h)*(delta_old(counter+1)delta_old(counter1));
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(N1))^2+...
b*sin(phi(N)delta_old(N))...
2*sqrt(ur*P)/(2*h)*(delta_old(1)delta_old(N1));
delta_new=H*y;
error=[error norm(delta_newdelta_old)];
delta_old=delta_new;
end;
delta=delta_new;
error
