use newton methode deal with a non linear (equation system)

[whhdlx]

I cant get that followed picture that what i need. S1(i) without result be shown on the matlab.
may be newton methode and funktion F and jacobi JF should be separate writen on the matlab.
but how to wirte a funtion for the matrix

%%%Front-fixing Method
%%implicit schema
%%Parameter ersteelung
X0=1; Xmax=2; K=1; r=0.1; vol=0.2; T=1;y1=1;error=0.1;N1=100;

%%Gitter erstellung
N=50;
M=50;
k=T/(N+1); %Timestep;
h=(Xmax-1)/(M+1); %stockstep;
X=[X0:h:Xmax];
%% Boundary Condition
% final condition
for j=1:M+2
f1(j,N+1)=0;
end
%freiei randwert
S1(N+1)=K;
%
for i=1:N+1
f1(M+2,i)=0;
end
%%%Definition von f1(1,i)und f1(2,i)
for i=1:N
f1(1,i)=K-S1(i);
f1(2,i)=K-(1+h)*S1(i);
end;

%%Berechnung S1,f1(1,i),f1(2,i)
%%Rechnung a1,b1,c1(Parameters) für model
for j=2:M+1
a1(j)=1+k/(h^2)*vol^2*X(j)^2+r*k;
for i=1:N
b1(j)=(-k)/(2*h^2)*vol^2*X(j)^2+k/(2*h)*X(j)*(r-(S1(i+1)-S1(i))/(k*S1(i)));
c1(j)=(-k)/(2*h^2)*vol^2*X(j)^2-k/(2*h)*X(j)*(r-(S1(i+1)-S1(i))/(k*S1(i)));
end
end


%%Set up Matrix A
A=zeros(M,M);
for j=2:M-1
A(j-1,j-1)=c1(j);
A(j,j-1)=a1(j+1);
A(j+1,j-1)=b1(j+2);
end
A(M-1,M-1)=c1(M);
A(M,M-1)=a1(M+1);
%
disp('es ist A')
%%Set up y
y=zeros(M,1);
for j=2:M
y(j-1,1)=f1(j+1,i)
end
y(M,1)=S1(i)
%%%Set up G
G=zeros(M,1);
for i=1:N
for j=3:M
G(j,1)=f1(j+1,i+1)
end
G(1,1)=f1(2,i+1)-b1(2)*(K-S1(i))-a1(2)*(K-(1+h)*S1(i));
G(2,1)=f1(3,i+1)-b1(3)*(K-(1+h)*S1(i));
end
%%%%Set up Nichtlinear gleichung F
F=A*y-G
disp('es ist F')
%%%erstellung Jacobbi Matrix von F
%% set up ableitung von a1,b1,c1__a2 b2 c2
for j=2:M+1
a2(j)=0;
for i=1:N
b2(j)=X(j)*S1(i+1)/(2*h*S1(i)^2);
c2(j)=-X(j)*S1(i+1)/(2*h*S1(i)^2);
end
end
%Set up ableitung von A nach S1: A1
A1=zeros(M,M);
for j=2:M-1
A1(j-1,j-1)=c2(j);
A1(j+1,j-1)=b2(j+2);
end
A1(M-1,M-1)=c2(M);
%%%set up Ableitung von G
G1=zeros(M,1);
for i=1:N
G1(1,1)=-(b2(2)*(K-S1(i))-b1(2))+a1(2)*(1+h);
G1(2,1)=-K*b2(3)+(1+h)*(b1(3)+b2(3)*S1(i));
end
%%Jacobi Matrix Fs
F1=A1*y-G1;
Fs=[A,F1];
disp('es ist Fs')
disp(Fs)

%%%inverse von Jacobbi Fs
%Fs1=inv(Fs);
%%%NEUTON Verfahren
%%% einer gegebenen Funktion y
y=y1;q=1;
while abs(F)>error; q=q+1;
if abs (Fs)<error||q>N1;
disp('Ableitung gleich null oder Anzahl N der Interationen überschritten');
break;
end;
y=y-F/Fs;
end;


disp;
disp('S1');
disp(S1);
plot(X,f1,1))

 

[hacksaw]

Is your question "how to display the matrix S1 data?"

 

[whhdlx]

I would like to use Implicit methode to deal with none parabolic system and use Newton's method
ich want solve A(S1(i))*f1,i) = G(S1(i)),
followed F(f1,i),S1(i))=A(S1(i))*f1,i)-G(S1(i))=0;
to be y=(f1(3,i) ,......,f1(M+1,i),S1(i));
and then with iterative Method(Newton's method) y(k+1)=y(k)-F(y(k))/J(F(y(k));
** J(F(y(k))) is Jacobi Matrix.
I need to find f1(1,i),....f1(M+1,i) and S1(i)
In the result to find f1,1) und S1(i), that´s mean in point 0 all of f1,1) and Value of S1(i) for all point of time.
And the diagram should the above showed picture.

 

[IRstuff]

Is this for school?

TTFN

FAQ731-376
Chinese prisoner wins Nobel Peace Prize

 

[hacksaw]

clearly it is not a home work assignment

as long the Jacobian is non-zero it should werk

 

[hacksaw]

How do you want this line of code to be evaluated?

a1(j)=1+k/(h^2)*vol^2*X(j)^2+r*k;

is it a1(j)=1+(k/(h^2))*vol^2*X(j)^2+r*k;

or a1(j)=1+k/(h^2*vol^2*X(j)^2)+r*k;

 

[hacksaw]

one additional comment:

"iterative Method(Newton's method) y(k+1)=y(k)-F(y(k))/J(F(y(k));"
where "J(F(y(k)))is Jacobi Matrix" is not correct,

it should be:

y(k+1)=y(k)-F(y(k))/det(J(F(y(k)))

 

[FeX32]

Why are you programming your own Newton-Euler method to solve a DE when you can program ode45 to do it with Runge-Kutta, which is better. The only time I programmed this manually was in an undergraduate numerical math class. hmmm.......


Fe

 

[hacksaw]

on closer examination, agreed